require(ggthemes)
require(ggplot2)
require(vegan)
require(dplyr)
require(data.table)
require(tidyr)

Reading in the new bird data

From Tyson’s csv. Creates a column for the abundance of birds by route (data is by stop).

setwd("/Users/Liv/Documents/NCEAS_GIT/NCEAS-RENCI_2014/BBS_data/")
new_bbs <- read.csv("fifty7_withRTENO.csv")
new_bbs$sum_route_abundance <- rowSums(new_bbs[grep("^Stop[0-9]+", names(new_bbs))])
head(new_bbs)
##   RouteDataID countrynum statenum Route RPID year AOU Stop1 Stop2 Stop3
## 1     6243648        124       62    16  101 2001  70     0     0     0
## 2     6243648        124       62    16  101 2001 100     0     0     0
## 3     6243648        124       62    16  101 2001 370     0     0     0
## 4     6243648        124       62    16  101 2001 380     0     0     0
## 5     6243648        124       62    16  101 2001 510     0     5     6
## 6     6243648        124       62    16  101 2001 710    46     0    19
##   Stop4 Stop5 Stop6 Stop7 Stop8 Stop9 Stop10 Stop11 Stop12 Stop13 Stop14
## 1     0     0     0     0     0     0      0      0      0      0      0
## 2     0     0     0     0     1     0      0      0      0      1      0
## 3     0     0     0     0     0     0      0      1      0      0      0
## 4     0     0     0     0     0     0      0      0      0      0      0
## 5     0     1     1    10     3     1      5      4      0      0      0
## 6    18     0     0     0     0     0      0      0      0      0      0
##   Stop15 Stop16 Stop17 Stop18 Stop19 Stop20 Stop21 Stop22 Stop23 Stop24
## 1      0      0      0      0      0      0      0      2      0      0
## 2      0      2      0      2      0      0      0      0      0      0
## 3      0      0      0      0      0      0      0      0      0      0
## 4      0      0      0      0      0      0      0      0      0      0
## 5      0      0      3      4      0      0      4      0      1      0
## 6      0      0      0      0      0      0      0      0      0      0
##   Stop25 Stop26 Stop27 Stop28 Stop29 Stop30 Stop31 Stop32 Stop33 Stop34
## 1      0      0      0      0      0      0      0      0      0      0
## 2      0      0      0      3      3      0      0      0      0      0
## 3      0      0      0      0      0      0      0      0      0      0
## 4      0      0      0      0      0      0      0      0      0      0
## 5      4      0      0      1      5      2      0      0      0      0
## 6      0      0      0      0      0      0      0      0      0      0
##   Stop35 Stop36 Stop37 Stop38 Stop39 Stop40 Stop41 Stop42 Stop43 Stop44
## 1      0      0      0      0      0      0      0      0      0      0
## 2      0      0      0      0      0      1      0      0      0      0
## 3      0      0      0      0      0      0      0      0      0      0
## 4      0      0      0      0      0      0      0      0      0      0
## 5      2      0      0      0      0      0      0      0      1      1
## 6      0      0      0      0      0      0      0      0      0      0
##   Stop45 Stop46 Stop47 Stop48 Stop49 Stop50 RTENO sum_route_abundance
## 1      0      0      0      0      0      1 62016                   3
## 2      0      0      0      0      0      0 62016                  13
## 3      0      0      0      0      0      0 62016                   1
## 4      1      0      0      0      0      0 62016                   1
## 5      0      0      0      0      0      1 62016                  65
## 6      0      0      0      0      0      0 62016                  83
names(new_bbs)
##  [1] "RouteDataID"         "countrynum"          "statenum"           
##  [4] "Route"               "RPID"                "year"               
##  [7] "AOU"                 "Stop1"               "Stop2"              
## [10] "Stop3"               "Stop4"               "Stop5"              
## [13] "Stop6"               "Stop7"               "Stop8"              
## [16] "Stop9"               "Stop10"              "Stop11"             
## [19] "Stop12"              "Stop13"              "Stop14"             
## [22] "Stop15"              "Stop16"              "Stop17"             
## [25] "Stop18"              "Stop19"              "Stop20"             
## [28] "Stop21"              "Stop22"              "Stop23"             
## [31] "Stop24"              "Stop25"              "Stop26"             
## [34] "Stop27"              "Stop28"              "Stop29"             
## [37] "Stop30"              "Stop31"              "Stop32"             
## [40] "Stop33"              "Stop34"              "Stop35"             
## [43] "Stop36"              "Stop37"              "Stop38"             
## [46] "Stop39"              "Stop40"              "Stop41"             
## [49] "Stop42"              "Stop43"              "Stop44"             
## [52] "Stop45"              "Stop46"              "Stop47"             
## [55] "Stop48"              "Stop49"              "Stop50"             
## [58] "RTENO"               "sum_route_abundance"
new_bbs <- new_bbs[ , c(1:7, 58, 59)] #selects everything except the individual stop abundances

Bird trait data

Read in and pare down the bird trait data to the scientific, common names for species and their habitat, nesting, diest, behaviour and conservation status, and 4&6 letter codes.

bird_traits <- read.csv("OH_BBS_1966-2013_traits.csv")
bird_traits <- bird_traits[ , c(1,2,3,4,5,6,7,8,9)]
head(bird_traits[2])
##         Scientific.name
## 1        Junco hyemalis
## 2      Colaptes auratus
## 3   Empidonax virescens
## 4     Empidonax alnorum
## 5 Botaurus lentiginosus
## 6         Anas rubripes
names(bird_traits) <- toupper(names(bird_traits))
names(bird_traits) <- gsub("[.]", "_", names(bird_traits))  # using the [] to get the period to be used literally. Or can use "fixed = TRUE"
names(bird_traits)
## [1] "X4_LETTER_CODE"  "SCIENTIFIC_NAME" "X6_LETTER_CODE"  "HABITAT"        
## [5] "DIET"            "NESTING"         "BEHAVIOR"        "CONSERVATION"   
## [9] "COMMON_NAME"
bird_traits$SCIENTIFIC_NAME <- toupper(bird_traits$SCIENTIFIC_NAME)
head(bird_traits)
##   X4_LETTER_CODE       SCIENTIFIC_NAME X6_LETTER_CODE       HABITAT
## 1           DEJU        JUNCO HYEMALIS         JUNHYE        forest
## 2           NOFL      COLAPTES AURATUS         COLAUT open_woodland
## 3           ACFL   EMPIDONAX VIRESCENS         EMPVIR        forest
## 4         ALFL       EMPIDONAX ALNORUM         EMPALN         scrub
## 5           AMBI BOTAURUS LENTIGINOSUS        BOTLEN          marsh
## 6          ABDU          ANAS RUBRIPES         ANARUB     lake_pond
##      DIET NESTING       BEHAVIOR CONSERVATION
## 1   seeds  ground ground_forager           LC
## 2 insects  cavity ground_forager           LC
## 3 insects    tree   fly_catching           LC
## 4 insects   shrub   fly_catching           LC
## 5    fish  ground       stalking           LC
## 6 insects  ground        dabbler           LC
##                                 COMMON_NAME
## 1     (Slate-colored Junco) Dark-eyed Junco
## 2 (Yellow-shafted Flicker) Northern Flicker
## 3                        Acadian Flycatcher
## 4                          Alder Flycatcher
## 5                          American Bittern
## 6                       American Black Duck

Bird code AOU data

BBS data comes with a different kind of AOU code. This chunk tidies up the AOU code dataframe and makes the common and scienitific names compatible with the BBS data. Need to use the common names, scientific names are slightly off in some cases.

AOU_codes <- read.csv("raw_data/AOU_codes.csv") ## a cleaned up (ie header removed, split into columns - no substantive changes to any species names or column headers) version of: ftp://ftpext.usgs.gov/pub/er/md/laurel/BBS/DataFiles/SpeciesList.txt

head(AOU_codes)
##   Seq   AOU          English_Common_Name           French_Common_Name
## 1   4 10010                Great Tinamou                Grand Tinamou
## 2   7 10030               Little Tinamou                 Tinamou soui
## 3  10 40080              Thicket Tinamou             Tinamou cannelle
## 4  13 10050       Slaty-breasted Tinamou           Tinamou de Boucard
## 5  16  1770 Black-bellied Whistling-Duck Dendrocygne \xe0 ventre noir
## 6  19 10200   West Indian Whistling-Duck     Dendrocygne des Antilles
##        Spanish_Common_Name        ORDER    Family        Genus     Species
## 1            Tinamus major Tinamiformes Tinamidae      Tinamus       major
## 2        Crypturellus soui Tinamiformes Tinamidae Crypturellus        soui
## 3 Crypturellus cinnamomeus Tinamiformes Tinamidae Crypturellus cinnamomeus
## 4    Crypturellus boucardi Tinamiformes Tinamidae Crypturellus    boucardi
## 5   Dendrocygna autumnalis Anseriformes  Anatidae  Dendrocygna  autumnalis
## 6      Dendrocygna arborea Anseriformes  Anatidae  Dendrocygna     arborea
AOU_codes$common_name <- toupper(AOU_codes$English_Common_Name)
names(AOU_codes)
##  [1] "Seq"                 "AOU"                 "English_Common_Name"
##  [4] "French_Common_Name"  "Spanish_Common_Name" "ORDER"              
##  [7] "Family"              "Genus"               "Species"            
## [10] "common_name"
AOU_codes <- AOU_codes[ , c(2, 10, 6:9)]
names(AOU_codes) <- toupper(names(AOU_codes))
head(AOU_codes)
##     AOU                  COMMON_NAME        ORDER    FAMILY        GENUS
## 1 10010                GREAT TINAMOU Tinamiformes Tinamidae      Tinamus
## 2 10030               LITTLE TINAMOU Tinamiformes Tinamidae Crypturellus
## 3 40080              THICKET TINAMOU Tinamiformes Tinamidae Crypturellus
## 4 10050       SLATY-BREASTED TINAMOU Tinamiformes Tinamidae Crypturellus
## 5  1770 BLACK-BELLIED WHISTLING-DUCK Anseriformes  Anatidae  Dendrocygna
## 6 10200   WEST INDIAN WHISTLING-DUCK Anseriformes  Anatidae  Dendrocygna
##       SPECIES
## 1       major
## 2        soui
## 3 cinnamomeus
## 4    boucardi
## 5  autumnalis
## 6     arborea
AOU_codes$SCIENTIFIC_NAME <- paste(AOU_codes$GENUS, AOU_codes$SPECIES, sep = " ")
AOU_codes$SCIENTIFIC_NAME <- toupper(AOU_codes$SCIENTIFIC_NAME)
#testing that all our traits at leat are in the new AOU codes:
bird_traits$SCIENTIFIC_NAME %in% AOU_codes$SCIENTIFIC_NAME # = not all true.
##   [1]  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE
##  [12]  TRUE  TRUE FALSE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE
##  [23]  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE FALSE  TRUE  TRUE  TRUE FALSE
##  [34]  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE
##  [45]  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE FALSE
##  [56]  TRUE  TRUE FALSE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE
##  [67]  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE FALSE
##  [78]  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE
##  [89]  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE
## [100]  TRUE FALSE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE
## [111]  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE
## [122]  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE
## [133]  TRUE  TRUE FALSE  TRUE  TRUE  TRUE FALSE FALSE  TRUE  TRUE  TRUE
## [144]  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE
## [155]  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE
## [166] FALSE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE FALSE FALSE FALSE FALSE
## [177] FALSE FALSE FALSE  TRUE  TRUE  TRUE  TRUE FALSE  TRUE  TRUE  TRUE
## [188]  TRUE FALSE  TRUE  TRUE  TRUE  TRUE FALSE  TRUE FALSE  TRUE  TRUE
## [199]  TRUE  TRUE FALSE  TRUE  TRUE  TRUE  TRUE
AOU_codes$SCIENTIFIC_NAME %in% bird_traits$SCIENTIFIC_NAME # many many falses
##    [1] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
##   [12] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
##   [23]  TRUE  TRUE FALSE  TRUE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
##   [34] FALSE FALSE FALSE  TRUE  TRUE FALSE FALSE FALSE FALSE  TRUE FALSE
##   [45] FALSE FALSE  TRUE FALSE FALSE FALSE FALSE FALSE FALSE  TRUE FALSE
##   [56] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
##   [67] FALSE FALSE FALSE FALSE FALSE FALSE FALSE  TRUE FALSE FALSE FALSE
##   [78]  TRUE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
##   [89] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE  TRUE FALSE FALSE
##  [100] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE  TRUE FALSE  TRUE
##  [111] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
##  [122] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
##  [133]  TRUE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
##  [144] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
##  [155] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
##  [166] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
##  [177] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
##  [188] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
##  [199] FALSE FALSE FALSE FALSE FALSE FALSE  TRUE FALSE FALSE FALSE FALSE
##  [210] FALSE FALSE FALSE FALSE  TRUE FALSE  TRUE FALSE  TRUE FALSE FALSE
##  [221]  TRUE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE  TRUE  TRUE
##  [232] FALSE FALSE  TRUE  TRUE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
##  [243] FALSE FALSE  TRUE  TRUE FALSE FALSE FALSE FALSE  TRUE FALSE FALSE
##  [254] FALSE FALSE FALSE FALSE FALSE FALSE  TRUE FALSE FALSE FALSE  TRUE
##  [265]  TRUE  TRUE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
##  [276] FALSE  TRUE  TRUE FALSE FALSE FALSE FALSE FALSE  TRUE FALSE FALSE
##  [287] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
##  [298] FALSE FALSE FALSE  TRUE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
##  [309] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE  TRUE
##  [320] FALSE FALSE FALSE  TRUE FALSE FALSE FALSE FALSE FALSE  TRUE FALSE
##  [331]  TRUE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
##  [342] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE  TRUE
##  [353] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE  TRUE FALSE
##  [364] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
##  [375] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
##  [386] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE  TRUE
##  [397] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
##  [408] FALSE FALSE FALSE  TRUE  TRUE FALSE FALSE FALSE FALSE FALSE FALSE
##  [419] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
##  [430] FALSE  TRUE FALSE FALSE FALSE  TRUE FALSE FALSE FALSE FALSE FALSE
##  [441] FALSE FALSE  TRUE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
##  [452] FALSE FALSE FALSE  TRUE  TRUE FALSE  TRUE FALSE FALSE FALSE FALSE
##  [463] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
##  [474] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
##  [485] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE  TRUE FALSE FALSE
##  [496] FALSE FALSE FALSE FALSE FALSE FALSE FALSE  TRUE FALSE FALSE FALSE
##  [507]  TRUE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
##  [518] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
##  [529] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
##  [540] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
##  [551] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
##  [562] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE  TRUE
##  [573] FALSE  TRUE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
##  [584] FALSE FALSE FALSE FALSE  TRUE FALSE FALSE FALSE FALSE FALSE FALSE
##  [595] FALSE FALSE FALSE  TRUE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
##  [606] FALSE FALSE FALSE FALSE FALSE  TRUE FALSE FALSE FALSE FALSE  TRUE
##  [617] FALSE FALSE FALSE FALSE FALSE FALSE  TRUE FALSE FALSE FALSE FALSE
##  [628] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
##  [639] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE  TRUE FALSE
##  [650] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
##  [661] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
##  [672] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
##  [683] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
##  [694] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
##  [705] FALSE FALSE FALSE FALSE FALSE  TRUE FALSE FALSE FALSE FALSE FALSE
##  [716] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
##  [727] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
##  [738] FALSE FALSE FALSE  TRUE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
##  [749] FALSE FALSE FALSE FALSE FALSE  TRUE FALSE FALSE FALSE FALSE FALSE
##  [760] FALSE FALSE  TRUE FALSE  TRUE FALSE FALSE FALSE FALSE FALSE FALSE
##  [771] FALSE  TRUE  TRUE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
##  [782] FALSE FALSE FALSE  TRUE FALSE FALSE FALSE FALSE  TRUE FALSE FALSE
##  [793] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
##  [804] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
##  [815] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
##  [826] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
##  [837] FALSE FALSE FALSE FALSE FALSE FALSE  TRUE FALSE FALSE  TRUE FALSE
##  [848] FALSE FALSE  TRUE  TRUE  TRUE FALSE FALSE  TRUE FALSE FALSE FALSE
##  [859] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE  TRUE FALSE
##  [870] FALSE FALSE FALSE FALSE FALSE FALSE FALSE  TRUE FALSE FALSE FALSE
##  [881] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
##  [892] FALSE FALSE FALSE  TRUE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
##  [903] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
##  [914] FALSE  TRUE FALSE FALSE  TRUE FALSE FALSE FALSE FALSE  TRUE FALSE
##  [925] FALSE FALSE  TRUE FALSE FALSE  TRUE FALSE FALSE FALSE  TRUE FALSE
##  [936]  TRUE  TRUE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
##  [947] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
##  [958] FALSE FALSE FALSE  TRUE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
##  [969] FALSE FALSE FALSE  TRUE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
##  [980]  TRUE FALSE FALSE  TRUE  TRUE FALSE FALSE FALSE FALSE FALSE FALSE
##  [991]  TRUE FALSE FALSE FALSE FALSE FALSE  TRUE FALSE  TRUE  TRUE FALSE
## [1002] FALSE  TRUE FALSE  TRUE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [1013] FALSE FALSE FALSE FALSE  TRUE FALSE FALSE FALSE FALSE  TRUE  TRUE
## [1024] FALSE FALSE  TRUE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [1035] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE  TRUE FALSE FALSE
## [1046]  TRUE  TRUE FALSE FALSE FALSE FALSE  TRUE FALSE  TRUE  TRUE FALSE
## [1057] FALSE FALSE FALSE FALSE  TRUE FALSE FALSE FALSE FALSE FALSE FALSE
## [1068] FALSE FALSE  TRUE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [1079] FALSE FALSE FALSE FALSE FALSE FALSE  TRUE FALSE FALSE FALSE FALSE
## [1090] FALSE FALSE FALSE FALSE FALSE FALSE  TRUE FALSE FALSE FALSE  TRUE
## [1101]  TRUE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [1112] FALSE  TRUE FALSE FALSE FALSE  TRUE FALSE  TRUE FALSE FALSE FALSE
## [1123] FALSE  TRUE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [1134] FALSE FALSE FALSE  TRUE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [1145] FALSE FALSE FALSE FALSE FALSE FALSE FALSE  TRUE FALSE FALSE FALSE
## [1156] FALSE FALSE FALSE FALSE FALSE FALSE FALSE  TRUE FALSE  TRUE  TRUE
## [1167]  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE FALSE FALSE FALSE FALSE FALSE
## [1178] FALSE FALSE FALSE FALSE FALSE FALSE  TRUE  TRUE FALSE FALSE FALSE
## [1189]  TRUE FALSE FALSE  TRUE  TRUE FALSE FALSE  TRUE  TRUE FALSE  TRUE
## [1200] FALSE  TRUE  TRUE FALSE FALSE FALSE FALSE  TRUE FALSE FALSE FALSE
## [1211]  TRUE  TRUE FALSE FALSE FALSE FALSE FALSE FALSE FALSE  TRUE FALSE
## [1222] FALSE FALSE FALSE  TRUE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [1233] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [1244] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [1255] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [1266] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE  TRUE FALSE
## [1277] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [1288] FALSE FALSE FALSE FALSE FALSE FALSE  TRUE FALSE  TRUE  TRUE FALSE
## [1299] FALSE  TRUE FALSE FALSE  TRUE  TRUE FALSE FALSE FALSE FALSE  TRUE
## [1310]  TRUE FALSE  TRUE FALSE FALSE FALSE FALSE FALSE FALSE FALSE  TRUE
## [1321] FALSE  TRUE FALSE  TRUE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [1332] FALSE FALSE  TRUE FALSE FALSE FALSE FALSE FALSE FALSE  TRUE FALSE
## [1343] FALSE FALSE FALSE FALSE FALSE FALSE FALSE  TRUE FALSE FALSE FALSE
## [1354]  TRUE FALSE FALSE FALSE FALSE FALSE FALSE  TRUE FALSE  TRUE FALSE
## [1365] FALSE FALSE FALSE FALSE  TRUE  TRUE  TRUE FALSE FALSE FALSE  TRUE
## [1376] FALSE FALSE  TRUE FALSE FALSE FALSE  TRUE FALSE FALSE FALSE FALSE
## [1387] FALSE FALSE  TRUE FALSE FALSE FALSE FALSE FALSE FALSE  TRUE FALSE
## [1398] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE  TRUE FALSE
## [1409] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [1420] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [1431] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [1442] FALSE FALSE FALSE FALSE  TRUE FALSE FALSE FALSE FALSE FALSE  TRUE
## [1453] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [1464] FALSE
head(bird_traits)
##   X4_LETTER_CODE       SCIENTIFIC_NAME X6_LETTER_CODE       HABITAT
## 1           DEJU        JUNCO HYEMALIS         JUNHYE        forest
## 2           NOFL      COLAPTES AURATUS         COLAUT open_woodland
## 3           ACFL   EMPIDONAX VIRESCENS         EMPVIR        forest
## 4         ALFL       EMPIDONAX ALNORUM         EMPALN         scrub
## 5           AMBI BOTAURUS LENTIGINOSUS        BOTLEN          marsh
## 6          ABDU          ANAS RUBRIPES         ANARUB     lake_pond
##      DIET NESTING       BEHAVIOR CONSERVATION
## 1   seeds  ground ground_forager           LC
## 2 insects  cavity ground_forager           LC
## 3 insects    tree   fly_catching           LC
## 4 insects   shrub   fly_catching           LC
## 5    fish  ground       stalking           LC
## 6 insects  ground        dabbler           LC
##                                 COMMON_NAME
## 1     (Slate-colored Junco) Dark-eyed Junco
## 2 (Yellow-shafted Flicker) Northern Flicker
## 3                        Acadian Flycatcher
## 4                          Alder Flycatcher
## 5                          American Bittern
## 6                       American Black Duck
bird_traits$COMMON_NAME %in% AOU_codes$English_Common_Name # = all true!
##   [1] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
##  [12] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
##  [23] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
##  [34] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
##  [45] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
##  [56] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
##  [67] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
##  [78] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
##  [89] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [100] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [111] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [122] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [133] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [144] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [155] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [166] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [177] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [188] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [199] FALSE FALSE FALSE FALSE FALSE FALSE FALSE
AOU_codes$English_Common_Name %in% bird_traits$COMMON_NAME # many many falses
## logical(0)
unique(new_bbs$AOU) %in% unique(AOU_codes$AOU) # good.
##   [1] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
##  [15] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
##  [29] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
##  [43] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
##  [57] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
##  [71] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
##  [85] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
##  [99] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [113] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [127] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [141] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [155] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [169] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [183] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [197] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [211] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [225] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [239] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [253] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [267] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [281] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [295] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [309] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [323] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [337] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [351] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [365] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [379] TRUE TRUE TRUE TRUE
# merge the AOU codes into the ohio data

Merging AOU codes and traits

Putting the trait and AOU codes together -

head(AOU_codes)
##     AOU                  COMMON_NAME        ORDER    FAMILY        GENUS
## 1 10010                GREAT TINAMOU Tinamiformes Tinamidae      Tinamus
## 2 10030               LITTLE TINAMOU Tinamiformes Tinamidae Crypturellus
## 3 40080              THICKET TINAMOU Tinamiformes Tinamidae Crypturellus
## 4 10050       SLATY-BREASTED TINAMOU Tinamiformes Tinamidae Crypturellus
## 5  1770 BLACK-BELLIED WHISTLING-DUCK Anseriformes  Anatidae  Dendrocygna
## 6 10200   WEST INDIAN WHISTLING-DUCK Anseriformes  Anatidae  Dendrocygna
##       SPECIES          SCIENTIFIC_NAME
## 1       major            TINAMUS MAJOR
## 2        soui        CRYPTURELLUS SOUI
## 3 cinnamomeus CRYPTURELLUS CINNAMOMEUS
## 4    boucardi    CRYPTURELLUS BOUCARDI
## 5  autumnalis   DENDROCYGNA AUTUMNALIS
## 6     arborea      DENDROCYGNA ARBOREA
head(new_bbs)
##   RouteDataID countrynum statenum Route RPID year AOU RTENO
## 1     6243648        124       62    16  101 2001  70 62016
## 2     6243648        124       62    16  101 2001 100 62016
## 3     6243648        124       62    16  101 2001 370 62016
## 4     6243648        124       62    16  101 2001 380 62016
## 5     6243648        124       62    16  101 2001 510 62016
## 6     6243648        124       62    16  101 2001 710 62016
##   sum_route_abundance
## 1                   3
## 2                  13
## 3                   1
## 4                   1
## 5                  65
## 6                  83
AOU_codes$COMMON_NAME <- AOU_codes$English_Common_Name

traits_aou <- inner_join(bird_traits, AOU_codes)
## Joining by: "SCIENTIFIC_NAME"
BBS_traits <- inner_join(new_bbs, traits_aou)
## Joining by: "AOU"
names(BBS_traits)
##  [1] "AOU"                 "RouteDataID"         "countrynum"         
##  [4] "statenum"            "Route"               "RPID"               
##  [7] "year"                "RTENO"               "sum_route_abundance"
## [10] "SCIENTIFIC_NAME"     "X4_LETTER_CODE"      "X6_LETTER_CODE"     
## [13] "HABITAT"             "DIET"                "NESTING"            
## [16] "BEHAVIOR"            "CONSERVATION"        "COMMON_NAME"        
## [19] "ORDER"               "FAMILY"              "GENUS"              
## [22] "SPECIES"

Neonics data by route

Based on Tyson’s GIS analysis.

neonics_buffers <- read.csv("../Pesticides/pest_buff_overtime.csv")
head(neonics_buffers)
##   YEAR RTENO buffer   COMPOUND high_kg_buff low_kg_buff ag_pix_buff
## 1 1992 66085   1000   ATRAZINE       859.34      847.89       52707
## 2 1992 66085   1000 GLYPHOSATE       108.95       77.76       52707
## 3 1992 66085  10000   ATRAZINE      8348.23     8223.01      502193
## 4 1992 66085  10000 GLYPHOSATE      1075.86      795.61      502193
## 5 1992 66085    200   ATRAZINE       181.23      178.84       11132
## 6 1992 66085    200 GLYPHOSATE        22.95       16.33       11132

Merge BBS_traits and neonics

sort(unique(neonics_buffers$YEAR))
##  [1] 1992 1993 1994 1995 1996 1997 1998 1999 2000 2001 2002 2003 2004 2005
## [15] 2006 2007 2008 2009 2010 2011
names(BBS_traits) <- toupper(names(BBS_traits))
sort(unique(BBS_traits$YEAR))
##  [1] 1995 1996 1997 1998 1999 2000 2001 2002 2003 2004 2005 2006 2007 2008
## [15] 2009 2010 2011 2012 2013

Therefore, we need to select the period 1995 - 2011 for merging.

BBS_neonics <- inner_join(BBS_traits[ BBS_traits$YEAR < 2012, ], neonics_buffers[ neonics_buffers$YEAR > 1994, ]) # joins by RTENO and YEAR which works. 
## Joining by: c("YEAR", "RTENO")

Landuse data by route

Based on Tyson’s GIS analysis. Notes for the csv: reclass is Anderson Level I land cover: (0=nodata, 1=water, 2=developed, 3=barren, 4=forest, 5=grassland/shrub, 6=agriculture, 7=wetlands). total is the number of pixels (multiply by 0.009 to get square km).

landuse_buffers <- read.csv("../Pesticides/LC_buffers_overtime.csv")
head(landuse_buffers)
##   BBS_route reclass buffer YEAR total_pix     km2 RTENO
## 1      2310       1   1000 1992         7  0.0063  2310
## 2      2310       2   1000 1992      7447  6.7023  2310
## 3      2310       3   1000 1992       106  0.0954  2310
## 4      2310       6   1000 1992     70656 63.5904  2310
## 5      2310       7   1000 1992      1155  1.0395  2310
## 6      2310       1  10000 1992      2435  2.1915  2310
reclass_table <- data.frame("reclass" = 0:7, "land_cover" = c("nodata", "water", "developed", "barren", "forest", "grassland_shrub", "agriculture", "wetlands"))

landuse_buffers <- merge(landuse_buffers, reclass_table, by = "reclass")
class(landuse_buffers$land_cover)
## [1] "factor"
head(landuse_buffers)
##   reclass BBS_route buffer YEAR total_pix    km2 RTENO land_cover
## 1       1      2310   1000 1992       7.0 0.0063  2310      water
## 2       1      2338    400 1996     229.8 0.2068 66049      water
## 3       1      2322   1000 1993     205.1 0.1846 66022      water
## 4       1      2359   2000 1992    1826.0 1.6434 66113      water
## 5       1      2371    400 1994     210.9 0.1898 66188      water
## 6       1      2310  10000 1992    2435.0 2.1915  2310      water
tail(landuse_buffers)
##       reclass BBS_route buffer YEAR total_pix    km2 RTENO land_cover
## 56235       7      2317   2000 1993    4806.2 4.3256 66011   wetlands
## 56236       7      2334   2000 2006     593.0 0.5337 66042   wetlands
## 56237       7      2367    400 1996       0.0 0.0000 66165   wetlands
## 56238       7      2312    400 2000    1560.4 1.4044 66002   wetlands
## 56239       7      2337   2000 2000     460.9 0.4148 66047   wetlands
## 56240       7      2372   2000 2005     116.0 0.1044 66251   wetlands
names(landuse_buffers)
## [1] "reclass"    "BBS_route"  "buffer"     "YEAR"       "total_pix" 
## [6] "km2"        "RTENO"      "land_cover"
landuse_buffers2 <- landuse_buffers[c(2, 3, 4, 5, 7, 8)]
names(landuse_buffers2)
## [1] "BBS_route"  "buffer"     "YEAR"       "total_pix"  "RTENO"     
## [6] "land_cover"
class(landuse_buffers2$total_pix)
## [1] "numeric"
str(landuse_buffers2)
## 'data.frame':    56240 obs. of  6 variables:
##  $ BBS_route : int  2310 2338 2322 2359 2371 2310 2322 2346 2326 2371 ...
##  $ buffer    : int  1000 400 1000 2000 400 10000 10000 2000 10000 5000 ...
##  $ YEAR      : int  1992 1996 1993 1992 1994 1992 1993 2010 1995 1994 ...
##  $ total_pix : num  7 230 205 1826 211 ...
##  $ RTENO     : int  2310 66049 66022 66113 66188 2310 66022 66073 66027 66188 ...
##  $ land_cover: Factor w/ 8 levels "agriculture",..: 7 7 7 7 7 7 7 7 7 7 ...
head(landuse_buffers2)
##   BBS_route buffer YEAR total_pix RTENO land_cover
## 1      2310   1000 1992       7.0  2310      water
## 2      2338    400 1996     229.8 66049      water
## 3      2322   1000 1993     205.1 66022      water
## 4      2359   2000 1992    1826.0 66113      water
## 5      2371    400 1994     210.9 66188      water
## 6      2310  10000 1992    2435.0  2310      water
landuse_buffers_wide <- spread(landuse_buffers2, land_cover, total_pix)
tail(landuse_buffers_wide, 50)
##      BBS_route buffer YEAR RTENO agriculture barren developed forest
## 8351      2379   2000 2002 66907       17962     NA      5161  63955
## 8352      2379   2000 2003 66907       17968     NA      5161  63949
## 8353      2379   2000 2004 66907       17975     NA      5161  63942
## 8354      2379   2000 2005 66907       17981     NA      5161  63936
## 8355      2379   2000 2006 66907       17988     NA      5161  63929
## 8356      2379   2000 2007 66907       17990     NA      5161  63927
## 8357      2379   2000 2008 66907       17992     NA      5161  63925
## 8358      2379   2000 2009 66907       17993     NA      5161  63923
## 8359      2379   2000 2010 66907       17995     NA      5161  63921
## 8360      2379   2000 2011 66907       17997     NA      5161  63919
## 8361      2379   5000 1992 66907       78416   0.00     17636 137612
## 8362      2379   5000 1993 66907       78231   0.00     17635 137774
## 8363      2379   5000 1994 66907       78046   0.00     17634 137936
## 8364      2379   5000 1995 66907       77861   0.00     17634 138097
## 8365      2379   5000 1996 66907       77676   0.00     17633 138259
## 8366      2379   5000 1997 66907       77491   0.00     17632 138421
## 8367      2379   5000 1998 66907       77306   0.00     17631 138583
## 8368      2379   5000 1999 66907       77121   0.00     17631 138744
## 8369      2379   5000 2000 66907       76936   0.00     17630 138906
## 8370      2379   5000 2001 66907       76751  13.00     17629 139068
## 8371      2379   5000 2002 66907       76758  13.00     17629 139061
## 8372      2379   5000 2003 66907       76764  13.00     17629 139055
## 8373      2379   5000 2004 66907       76771  13.00     17629 139048
## 8374      2379   5000 2005 66907       76777  13.00     17629 139042
## 8375      2379   5000 2006 66907       76784  13.00     17629 139035
## 8376      2379   5000 2007 66907       76805  13.00     17636 139003
## 8377      2379   5000 2008 66907       76827  13.00     17642 138971
## 8378      2379   5000 2009 66907       76848  13.00     17649 138938
## 8379      2379   5000 2010 66907       76870  13.00     17655 138906
## 8380      2379   5000 2011 66907       76891  50.00     17662 138874
## 8381      2379  10000 1992 66907      240727  48.00     43000 326098
## 8382      2379  10000 1993 66907      240348  44.78     42997 326434
## 8383      2379  10000 1994 66907      239969  41.56     42994 326770
## 8384      2379  10000 1995 66907      239589  38.33     42990 327106
## 8385      2379  10000 1996 66907      239210  35.11     42987 327442
## 8386      2379  10000 1997 66907      238831  31.89     42984 327778
## 8387      2379  10000 1998 66907      238452  28.67     42981 328114
## 8388      2379  10000 1999 66907      238072  25.44     42977 328450
## 8389      2379  10000 2000 66907      237693  22.22     42974 328786
## 8390      2379  10000 2001 66907      237314  24.00     42971 329122
## 8391      2379  10000 2002 66907      237320  24.00     42971 329111
## 8392      2379  10000 2003 66907      237327  24.00     42971 329101
## 8393      2379  10000 2004 66907      237333  24.00     42971 329090
## 8394      2379  10000 2005 66907      237340  24.00     42971 329080
## 8395      2379  10000 2006 66907      237346  24.00     42971 329069
## 8396      2379  10000 2007 66907      237462  31.40     42987 328941
## 8397      2379  10000 2008 66907      237577  38.80     43003 328813
## 8398      2379  10000 2009 66907      237693  46.20     43018 328685
## 8399      2379  10000 2010 66907      237808  53.60     43034 328557
## 8400      2379  10000 2011 66907      237924  61.00     43050 328429
##      grassland_shrub water wetlands
## 8351             822  3285    144.0
## 8352             822  3285    144.0
## 8353             822  3285    144.0
## 8354             822  3285    144.0
## 8355             822  3285    144.0
## 8356             822  3285    144.2
## 8357             822  3285    144.4
## 8358             822  3285    144.6
## 8359             822  3285    144.8
## 8360             822  3285    145.0
## 8361            2307  4445    121.0
## 8362            2318  4457    121.8
## 8363            2329  4469    122.6
## 8364            2340  4482    123.3
## 8365            2351  4494    124.1
## 8366            2362  4506    124.9
## 8367            2373  4518    125.7
## 8368            2384  4531    126.4
## 8369            2395  4543    127.2
## 8370            2490  4578    611.0
## 8371            2490  4578    611.0
## 8372            2490  4578    611.0
## 8373            2490  4578    611.0
## 8374            2490  4578    611.0
## 8375            2490  4578    611.0
## 8376            2482  4580    613.6
## 8377            2474  4583    616.2
## 8378            2465  4585    618.8
## 8379            2457  4588    621.4
## 8380            2449  4590    624.0
## 8381            6283  6532    752.0
## 8382            6310  6550    756.4
## 8383            6338  6568    760.9
## 8384            6365  6585    765.3
## 8385            6393  6603    769.8
## 8386            6420  6621    774.2
## 8387            6448  6639    778.7
## 8388            6475  6656    783.1
## 8389            6503  6674    787.6
## 8390            6530  6692   1990.0
## 8391            6530  6700   1986.0
## 8392            6530  6708   1982.0
## 8393            6530  6717   1978.0
## 8394            6530  6725   1974.0
## 8395            6530  6733   1970.0
## 8396            6503  6742   1976.8
## 8397            6476  6752   1983.6
## 8398            6449  6761   1990.4
## 8399            6422  6771   1997.2
## 8400            6395  6780   2004.0
names(landuse_buffers_wide)
##  [1] "BBS_route"       "buffer"          "YEAR"           
##  [4] "RTENO"           "agriculture"     "barren"         
##  [7] "developed"       "forest"          "grassland_shrub"
## [10] "water"           "wetlands"

Merging BBS_neonics with landuse

sort(unique(landuse_buffers$YEAR)) # 1992 - 2011
##  [1] 1992 1993 1994 1995 1996 1997 1998 1999 2000 2001 2002 2003 2004 2005
## [15] 2006 2007 2008 2009 2010 2011
BBS_ln <- inner_join(BBS_neonics, landuse_buffers_wide[ landuse_buffers_wide$YEAR > 1994, ])
## Joining by: c("YEAR", "RTENO", "buffer")
# joins by: c("YEAR", "RTENO", "BBS_route", "buffer")

tail(BBS_ln)
##         YEAR RTENO buffer  AOU ROUTEDATAID COUNTRYNUM STATENUM ROUTE RPID
## 1429884 2002 66907   2000 7660     6245008        840       66   907  101
## 1429885 2002 66907    400 7660     6245008        840       66   907  101
## 1429886 2002 66907    400 7660     6245008        840       66   907  101
## 1429887 2002 66907    400 7660     6245008        840       66   907  101
## 1429888 2002 66907    400 7660     6245008        840       66   907  101
## 1429889 2002 66907    400 7660     6245008        840       66   907  101
##         SUM_ROUTE_ABUNDANCE SCIENTIFIC_NAME X4_LETTER_CODE X6_LETTER_CODE
## 1429884                   3   SIALIA SIALIS          EABL          SIASIA
## 1429885                   3   SIALIA SIALIS          EABL          SIASIA
## 1429886                   3   SIALIA SIALIS          EABL          SIASIA
## 1429887                   3   SIALIA SIALIS          EABL          SIASIA
## 1429888                   3   SIALIA SIALIS          EABL          SIASIA
## 1429889                   3   SIALIA SIALIS          EABL          SIASIA
##           HABITAT    DIET NESTING       BEHAVIOR CONSERVATION
## 1429884 grassland insects  cavity ground_forager           LC
## 1429885 grassland insects  cavity ground_forager           LC
## 1429886 grassland insects  cavity ground_forager           LC
## 1429887 grassland insects  cavity ground_forager           LC
## 1429888 grassland insects  cavity ground_forager           LC
## 1429889 grassland insects  cavity ground_forager           LC
##              COMMON_NAME         ORDER   FAMILY  GENUS SPECIES
## 1429884 Eastern Bluebird Passeriformes Turdidae Sialia  sialis
## 1429885 Eastern Bluebird Passeriformes Turdidae Sialia  sialis
## 1429886 Eastern Bluebird Passeriformes Turdidae Sialia  sialis
## 1429887 Eastern Bluebird Passeriformes Turdidae Sialia  sialis
## 1429888 Eastern Bluebird Passeriformes Turdidae Sialia  sialis
## 1429889 Eastern Bluebird Passeriformes Turdidae Sialia  sialis
##             COMPOUND high_kg_buff low_kg_buff ag_pix_buff BBS_route
## 1429884 THIAMETHOXAM    2.265e-04   2.265e-04        2186      2379
## 1429885     ATRAZINE    2.287e+01   2.269e+01        1470      2379
## 1429886     FIPRONIL    3.360e-02          NA        1470      2379
## 1429887   GLYPHOSATE    3.253e+01   3.199e+01        1470      2379
## 1429888 IMIDACLOPRID    2.421e-02   2.421e-02        1470      2379
## 1429889 THIAMETHOXAM    6.218e-07   6.218e-07           6      2379
##         agriculture barren developed forest grassland_shrub water wetlands
## 1429884       17962     NA      5161  63955             822  3285      144
## 1429885        1463     NA      1786  24349             300   740       17
## 1429886        1463     NA      1786  24349             300   740       17
## 1429887        1463     NA      1786  24349             300   740       17
## 1429888        1463     NA      1786  24349             300   740       17
## 1429889        1463     NA      1786  24349             300   740       17

Maybe a baby model?

names(BBS_ln)
##  [1] "YEAR"                "RTENO"               "buffer"             
##  [4] "AOU"                 "ROUTEDATAID"         "COUNTRYNUM"         
##  [7] "STATENUM"            "ROUTE"               "RPID"               
## [10] "SUM_ROUTE_ABUNDANCE" "SCIENTIFIC_NAME"     "X4_LETTER_CODE"     
## [13] "X6_LETTER_CODE"      "HABITAT"             "DIET"               
## [16] "NESTING"             "BEHAVIOR"            "CONSERVATION"       
## [19] "COMMON_NAME"         "ORDER"               "FAMILY"             
## [22] "GENUS"               "SPECIES"             "COMPOUND"           
## [25] "high_kg_buff"        "low_kg_buff"         "ag_pix_buff"        
## [28] "BBS_route"           "agriculture"         "barren"             
## [31] "developed"           "forest"              "grassland_shrub"    
## [34] "water"               "wetlands"
unique(BBS_ln$buffer)
## [1]  1000 10000   200  2000   400  5000
BBS_ln$agriculture_scaled <- scale(BBS_ln$agriculture)
BBS_ln$high_kg_buff_scaled <- scale(BBS_ln$high_kg_buff)
head(BBS_ln)
##   YEAR RTENO buffer  AOU ROUTEDATAID COUNTRYNUM STATENUM ROUTE RPID
## 1 1997 66001   1000 1320     6228264        840       66     1  101
## 2 1997 66001   1000 1320     6228264        840       66     1  101
## 3 1997 66001   1000 1320     6228264        840       66     1  101
## 4 1997 66001  10000 1320     6228264        840       66     1  101
## 5 1997 66001  10000 1320     6228264        840       66     1  101
## 6 1997 66001  10000 1320     6228264        840       66     1  101
##   SUM_ROUTE_ABUNDANCE    SCIENTIFIC_NAME X4_LETTER_CODE X6_LETTER_CODE
## 1                   2 ANAS PLATYRHYNCHOS           MALL         ANAPLA
## 2                   2 ANAS PLATYRHYNCHOS           MALL         ANAPLA
## 3                   2 ANAS PLATYRHYNCHOS           MALL         ANAPLA
## 4                   2 ANAS PLATYRHYNCHOS           MALL         ANAPLA
## 5                   2 ANAS PLATYRHYNCHOS           MALL         ANAPLA
## 6                   2 ANAS PLATYRHYNCHOS           MALL         ANAPLA
##     HABITAT  DIET NESTING BEHAVIOR CONSERVATION COMMON_NAME        ORDER
## 1 lake_pond seeds  ground  dabbler           LC     Mallard Anseriformes
## 2 lake_pond seeds  ground  dabbler           LC     Mallard Anseriformes
## 3 lake_pond seeds  ground  dabbler           LC     Mallard Anseriformes
## 4 lake_pond seeds  ground  dabbler           LC     Mallard Anseriformes
## 5 lake_pond seeds  ground  dabbler           LC     Mallard Anseriformes
## 6 lake_pond seeds  ground  dabbler           LC     Mallard Anseriformes
##     FAMILY GENUS       SPECIES     COMPOUND high_kg_buff low_kg_buff
## 1 Anatidae  Anas platyrhynchos     ATRAZINE           NA          NA
## 2 Anatidae  Anas platyrhynchos   GLYPHOSATE           NA          NA
## 3 Anatidae  Anas platyrhynchos IMIDACLOPRID           NA          NA
## 4 Anatidae  Anas platyrhynchos     ATRAZINE           NA          NA
## 5 Anatidae  Anas platyrhynchos   GLYPHOSATE           NA          NA
## 6 Anatidae  Anas platyrhynchos IMIDACLOPRID           NA          NA
##   ag_pix_buff BBS_route agriculture barren developed forest
## 1          NA      2311       75241      0      7387      0
## 2          NA      2311       75241      0      7387      0
## 3          NA      2311       75241      0      7387      0
## 4          NA      2311      997407   1057     96264  24625
## 5          NA      2311      997407   1057     96264  24625
## 6          NA      2311      997407   1057     96264  24625
##   grassland_shrub   water wetlands agriculture_scaled high_kg_buff_scaled
## 1               0   13.44        0            -0.3525                  NA
## 2               0   13.44        0            -0.3525                  NA
## 3               0   13.44        0            -0.3525                  NA
## 4            8417 8353.56     1241             3.5449                  NA
## 5            8417 8353.56     1241             3.5449                  NA
## 6            8417 8353.56     1241             3.5449                  NA
# test_mod2 <- glmer(data = BBS_ln[ BBS_ln$buffer == 200, ],
#                    family = poisson,
#                    formula= SUM_ROUTE_ABUNDANCE ~ agriculture_scaled + high_kg_buff_scaled:COMPOUND + DIET + (1|BBS_route) + (1|YEAR))
#summary(test_mod2)
require(coefplot2)
#coefplot2(test_mod2)


require(dismo)
## Loading required package: dismo
## Loading required package: raster
## Loading required package: sp
## 
## Attaching package: 'raster'
## 
## The following object is masked from 'package:tidyr':
## 
##     extract
## 
## The following object is masked from 'package:dplyr':
## 
##     select
BBS_200 <- BBS_ln[ BBS_ln$buffer == 200, ]
names(BBS_200)
##  [1] "YEAR"                "RTENO"               "buffer"             
##  [4] "AOU"                 "ROUTEDATAID"         "COUNTRYNUM"         
##  [7] "STATENUM"            "ROUTE"               "RPID"               
## [10] "SUM_ROUTE_ABUNDANCE" "SCIENTIFIC_NAME"     "X4_LETTER_CODE"     
## [13] "X6_LETTER_CODE"      "HABITAT"             "DIET"               
## [16] "NESTING"             "BEHAVIOR"            "CONSERVATION"       
## [19] "COMMON_NAME"         "ORDER"               "FAMILY"             
## [22] "GENUS"               "SPECIES"             "COMPOUND"           
## [25] "high_kg_buff"        "low_kg_buff"         "ag_pix_buff"        
## [28] "BBS_route"           "agriculture"         "barren"             
## [31] "developed"           "forest"              "grassland_shrub"    
## [34] "water"               "wetlands"            "agriculture_scaled" 
## [37] "high_kg_buff_scaled"
class(BBS_200$SPECIES)
## [1] "factor"
# gbm_neonic <- gbm.step(data = BBS_200, 
#                      gbm.x=c(1, 3, 15, 16, 17, 18, 23, 25, 26, 29:35), 
#                      gbm.y = 11, family = "poisson",
#                      tree.complexity = 14) #, 
#                       #step.size = 20, n.trees = 50,
#                      #learning.rate = 0.001, bag.fraction = 0.75)
# summary(gbm_neonic)
# 
# glmm_neonic <- glmer(data = BBS_200,
#                      formula = SUM_ROUTE_ABUNDANCE ~ high_kg_buff:COMPOUND + 
#                        high_kg_buff:forest:COMPOUND + 
#                        high_kg_buff:agriculture:COMPOUND + 
#                        forest + agriculture +
#                        (1|GENUS) + (1|YEAR),
#                      family = "poisson")
# 
# summary(glmm_neonic)
# coefplot2(glmm_neonic)
# 
# neonic_int<-gbm.interactions(gbm_neonic)
# neonic_int$interactions
# 
# gbm.perspec(gbm_neonic, 1,3)
# gbm.perspec(gbm_neonic, 2,3)
# gbm.perspec(gbm_neonic, 4,8)
# gbm.plot(gbm_neonic, 3,4)
# 
# save.image(file = "big_gbm.Rdata")
# load(file = "big_gbm.Rdata")
#saveRDS(gbm_neonic, file = "big_serialised_gbm.rds")

Some exploration!

head(landuse_buffers)
##   reclass BBS_route buffer YEAR total_pix    km2 RTENO land_cover
## 1       1      2310   1000 1992       7.0 0.0063  2310      water
## 2       1      2338    400 1996     229.8 0.2068 66049      water
## 3       1      2322   1000 1993     205.1 0.1846 66022      water
## 4       1      2359   2000 1992    1826.0 1.6434 66113      water
## 5       1      2371    400 1994     210.9 0.1898 66188      water
## 6       1      2310  10000 1992    2435.0 2.1915  2310      water
ggplot(landuse_buffers[landuse_buffers$land_cover %in% c("forest", "agriculture", "developed", "grassland_shrub") & landuse_buffers$buffer == 400, ], aes(x = YEAR, y = total_pix)) +
  geom_smooth(aes(colour = land_cover)) +
#  geom_density2d(aes(colour = land_cover)) +
  geom_point(aes(colour = land_cover), alpha = 0.3, position = position_jitter(width = 0.3))+
  theme_bw()+
  scale_colour_colorblind("Land Cover", labels = c("Agriculture", "Developed", "Forest", "Grass/Shrub"))+
  guides(colour=guide_legend(override.aes=list(fill=NA)))+
  theme(legend.position = "bottom", legend.key = element_blank(), legend.background = element_rect(colour = "grey"))
## geom_smooth: method="auto" and size of largest group is >=1000, so using gam with formula: y ~ s(x, bs = "cs"). Use 'method = x' to change the smoothing method.

plot of chunk unnamed-chunk-13

ggplot(landuse_buffers[landuse_buffers$land_cover %in% c("forest", "agriculture", "developed", "grassland_shrub") & landuse_buffers$buffer == 400, ], aes(x = YEAR, y = total_pix)) +
  geom_smooth(aes(colour = land_cover)) +
  geom_density2d(aes(colour = land_cover)) +
  #geom_point(aes(colour = land_cover), alpha = 0.3, position = position_jitter(width = 0.3))+
  theme_bw()+
  scale_colour_colorblind("Land Cover", labels = c("Agriculture", "Developed", "Forest", "Grass/Shrub"))+
  guides(colour = FALSE) +
  facet_wrap(~ land_cover, scales = "free")
## geom_smooth: method="auto" and size of largest group is >=1000, so using gam with formula: y ~ s(x, bs = "cs"). Use 'method = x' to change the smoothing method.
## geom_smooth: method="auto" and size of largest group is >=1000, so using gam with formula: y ~ s(x, bs = "cs"). Use 'method = x' to change the smoothing method.
## geom_smooth: method="auto" and size of largest group is >=1000, so using gam with formula: y ~ s(x, bs = "cs"). Use 'method = x' to change the smoothing method.
## geom_smooth: method="auto" and size of largest group is >=1000, so using gam with formula: y ~ s(x, bs = "cs"). Use 'method = x' to change the smoothing method.

plot of chunk unnamed-chunk-13

ggplot(landuse_buffers[landuse_buffers$land_cover %in% c("forest", "agriculture", "developed", "grassland_shrub") & landuse_buffers$buffer == 400 & landuse_buffers$YEAR %in% c(1995, 2000, 2005, 2010), ], aes(x = land_cover, y = total_pix)) +
 geom_violin(aes(fill = land_cover))+
  theme_bw()+
  scale_fill_colorblind("Land Cover", labels = c("Agriculture", "Developed", "Forest", "Grass/Shrub")) +
  facet_wrap(~YEAR)

plot of chunk unnamed-chunk-13

ggplot(landuse_buffers[landuse_buffers$land_cover %in% c("forest", "agriculture", "developed", "grassland_shrub") & landuse_buffers$buffer == 400 & landuse_buffers$YEAR %in% c(1995, 2000, 2005, 2010), ], aes(x = factor(YEAR), y = total_pix)) +
 geom_violin(aes(fill = land_cover))+
  theme_bw()+
  scale_fill_colorblind("Land Cover", labels = c("Agriculture", "Developed", "Forest", "Grass/Shrub")) +
  guides(fill = FALSE)+
  facet_wrap(~land_cover, scales = "free")

plot of chunk unnamed-chunk-13

summary_landuse <- landuse_buffers[landuse_buffers$land_cover %in% c("forest", "agriculture", "developed", "grassland_shrub"), ] %>% group_by(YEAR, land_cover, buffer) %>% filter(buffer == 400) %>% summarise(total_pix = sum(total_pix))

ggplot(summary_landuse, aes(x = YEAR, y = total_pix)) +
  geom_area(aes(fill = land_cover, colour = land_cover), position = "fill")+
  theme_bw()+
  scale_fill_manual("Land Cover", values = c("sienna4", "gray20", "forestgreen", "orange"), labels = c("Agriculture", "Developed", "Forest", "Grass/Shrub")) +
  scale_colour_manual("Land Cover", values = c("sienna4", "gray20", "forestgreen", "orange"), labels = c("Agriculture", "Developed", "Forest", "Grass/Shrub")) +
  theme(legend.position = "bottom", legend.background = element_rect(colour = "grey"))

plot of chunk unnamed-chunk-13

unique(landuse_buffers$buffer)
## [1]  1000   400  2000 10000  5000   200
summary_landuse3 <- landuse_buffers[landuse_buffers$land_cover %in% c("forest", "agriculture", "developed", "grassland_shrub"), ] %>% group_by(YEAR, land_cover, buffer) %>% summarise(total_pix = sum(total_pix))

ggplot(summary_landuse3, aes(x = YEAR, y = total_pix)) +
  geom_area(aes(fill = land_cover, colour = land_cover), position = "fill")+
  theme_bw()+
  scale_fill_manual("Land Cover", values = c("sienna4", "gray20", "forestgreen", "orange"), labels = c("Agriculture", "Developed", "Forest", "Grass/Shrub")) +
  scale_colour_manual("Land Cover", values = c("sienna4", "gray20", "forestgreen", "orange"), labels = c("Agriculture", "Developed", "Forest", "Grass/Shrub")) +
  theme(legend.position = "bottom", legend.background = element_rect(colour = "grey")) +
  facet_wrap(~buffer)

plot of chunk unnamed-chunk-13

# summary_landuse2 <- landuse_buffers%>% group_by(YEAR, land_cover, buffer) %>% filter(buffer == 400) %>% summarise(total_pix = sum(total_pix))
# 
# ggplot(summary_landuse2, aes(x = YEAR, y = total_pix)) +
#   geom_area(aes(fill = land_cover, colour = land_cover), position = "fill")+
#   theme_bw()+
#   #scale_fill_manual("Land Cover", values = c("sienna4", "gray20", "forestgreen", "orange"), labels = c("Agriculture", "Developed", "Forest", "Grass/Shrub")) +
#   #scale_colour_manual("Land Cover", values = c("sienna4", "gray20", "forestgreen", "orange"), labels = c("Agriculture", "Developed", "Forest", "Grass/Shrub")) +
#   theme(legend.position = "bottom", legend.background = element_rect(colour = "grey"))

Exploring the neonics

head(neonics_buffers)
##   YEAR RTENO buffer   COMPOUND high_kg_buff low_kg_buff ag_pix_buff
## 1 1992 66085   1000   ATRAZINE       859.34      847.89       52707
## 2 1992 66085   1000 GLYPHOSATE       108.95       77.76       52707
## 3 1992 66085  10000   ATRAZINE      8348.23     8223.01      502193
## 4 1992 66085  10000 GLYPHOSATE      1075.86      795.61      502193
## 5 1992 66085    200   ATRAZINE       181.23      178.84       11132
## 6 1992 66085    200 GLYPHOSATE        22.95       16.33       11132
head(neonics_buffers[is.na(neonics_buffers$high_kg_buff) | is.na(neonics_buffers$ag_pix_buff) | is.na(neonics_buffers$low_kg_buff), ])
##     YEAR RTENO buffer     COMPOUND high_kg_buff low_kg_buff ag_pix_buff
## 141 1995 66907  10000 IMIDACLOPRID       0.4244          NA      205449
## 144 1995 66907   5000 IMIDACLOPRID       0.1192          NA       70705
## 195 1997 66085   1000 IMIDACLOPRID       1.7907          NA       52402
## 198 1997 66085  10000 IMIDACLOPRID      16.5500          NA      499972
## 201 1997 66085    200 IMIDACLOPRID       0.3804          NA       11106
## 204 1997 66085   2000 IMIDACLOPRID       3.2648          NA       95762
head(neonics_buffers[is.na(neonics_buffers$high_kg_buff) | is.na(neonics_buffers$ag_pix_buff) , ])
##      YEAR RTENO buffer   COMPOUND high_kg_buff low_kg_buff ag_pix_buff
## 2135 1992 66011    200   ATRAZINE           NA          NA          NA
## 2136 1992 66011    200 GLYPHOSATE           NA          NA          NA
## 2149 1992 66908    200   ATRAZINE           NA          NA          NA
## 2150 1992 66908    200 GLYPHOSATE           NA          NA          NA
## 2167 1993 66011    200   ATRAZINE           NA          NA          NA
## 2168 1993 66011    200 GLYPHOSATE           NA          NA          NA
dim(neonics_buffers[is.na(neonics_buffers$high_kg_buff) | is.na(neonics_buffers$ag_pix_buff) | is.na(neonics_buffers$low_kg_buff), ])
## [1] 11997     7
dim(neonics_buffers[is.na(neonics_buffers$high_kg_buff) | is.na(neonics_buffers$ag_pix_buff), ])
## [1] 1501    7
# unique(neonics_buffers$YEAR)
# class(neonics_buffers$high_kg_buff)
# noyear_neonics_buffers <- neonics_buffers %>% 
#   group_by(BBS_route, buffer) %>% 
#   summarise(high_pest = sum(high_kg_buff), agpix = sum(ag_pix_buff))
# head(noyear_neonics_buffers)
# # spps <- spps[!is.na(spps$NMDS1) & !is.na(spps$NMDS2),]
# noyear_neonics_buffers <- noyear_neonics_buffers[!is.na(noyear_neonics_buffers$high_pest) & !is.na(noyear_neonics_buffers$agpix), ]
# 
# ggplot(noyear_neonics_buffers, aes(y = as.numeric(high_pest/agpix), x = as.numeric(buffer))) +
#   geom_line(aes(group = BBS_route)) +
#   geom_smooth()+
#   theme_bw()
# 
# ggplot(noyear_neonics_buffers, aes(y = high_pest, x = buffer)) +
#   geom_line(aes(group = BBS_route)) +
#   geom_smooth()+
#   theme_bw()
# 
# ggplot(noyear_neonics_buffers, aes(y = agpix, x = buffer)) +
#   geom_line(aes(group = BBS_route)) +
#   geom_smooth()+
#   theme_bw()
# 
# 
# # now to see if intensity has increased over tine
# names(neonics_buffers)
# year_neonics_buffers <- neonics_buffers %>%
#   group_by(BBS_route, YEAR, COMPOUND, buffer) %>%
#   summarise(intensity = high_kg_buff/ag_pix_buff) 
# 
# names(year_neonics_buffers)
# 
# ggplot(year_neonics_buffers, aes(y = intensity, x = YEAR)) +
#   geom_line(aes(group = BBS_route)) +
#   facet_grid(buffer~ COMPOUND, scale = "free") +
#   theme_bw()
# 
# year_neonics_buffers2 <- neonics_buffers %>%
#   group_by(BBS_route, YEAR, COMPOUND) %>%
#   summarise(intensity = mean(high_kg_buff/ag_pix_buff)) 
# 
# ggplot(year_neonics_buffers2, aes(y = intensity, x = YEAR)) + 
#   geom_line(aes(group = BBS_route)) +
#   facet_wrap(~ COMPOUND, scale = "free") +
#   theme_bw()
# 
# ggplot(year_neonics_buffers2, aes(y = intensity, x = YEAR)) + 
#   geom_line(aes(group = BBS_route)) +
#   facet_wrap(~ COMPOUND) +
#   theme_bw()

Using the big dataset: BBS_ln

head(BBS_ln)
##   YEAR RTENO buffer  AOU ROUTEDATAID COUNTRYNUM STATENUM ROUTE RPID
## 1 1997 66001   1000 1320     6228264        840       66     1  101
## 2 1997 66001   1000 1320     6228264        840       66     1  101
## 3 1997 66001   1000 1320     6228264        840       66     1  101
## 4 1997 66001  10000 1320     6228264        840       66     1  101
## 5 1997 66001  10000 1320     6228264        840       66     1  101
## 6 1997 66001  10000 1320     6228264        840       66     1  101
##   SUM_ROUTE_ABUNDANCE    SCIENTIFIC_NAME X4_LETTER_CODE X6_LETTER_CODE
## 1                   2 ANAS PLATYRHYNCHOS           MALL         ANAPLA
## 2                   2 ANAS PLATYRHYNCHOS           MALL         ANAPLA
## 3                   2 ANAS PLATYRHYNCHOS           MALL         ANAPLA
## 4                   2 ANAS PLATYRHYNCHOS           MALL         ANAPLA
## 5                   2 ANAS PLATYRHYNCHOS           MALL         ANAPLA
## 6                   2 ANAS PLATYRHYNCHOS           MALL         ANAPLA
##     HABITAT  DIET NESTING BEHAVIOR CONSERVATION COMMON_NAME        ORDER
## 1 lake_pond seeds  ground  dabbler           LC     Mallard Anseriformes
## 2 lake_pond seeds  ground  dabbler           LC     Mallard Anseriformes
## 3 lake_pond seeds  ground  dabbler           LC     Mallard Anseriformes
## 4 lake_pond seeds  ground  dabbler           LC     Mallard Anseriformes
## 5 lake_pond seeds  ground  dabbler           LC     Mallard Anseriformes
## 6 lake_pond seeds  ground  dabbler           LC     Mallard Anseriformes
##     FAMILY GENUS       SPECIES     COMPOUND high_kg_buff low_kg_buff
## 1 Anatidae  Anas platyrhynchos     ATRAZINE           NA          NA
## 2 Anatidae  Anas platyrhynchos   GLYPHOSATE           NA          NA
## 3 Anatidae  Anas platyrhynchos IMIDACLOPRID           NA          NA
## 4 Anatidae  Anas platyrhynchos     ATRAZINE           NA          NA
## 5 Anatidae  Anas platyrhynchos   GLYPHOSATE           NA          NA
## 6 Anatidae  Anas platyrhynchos IMIDACLOPRID           NA          NA
##   ag_pix_buff BBS_route agriculture barren developed forest
## 1          NA      2311       75241      0      7387      0
## 2          NA      2311       75241      0      7387      0
## 3          NA      2311       75241      0      7387      0
## 4          NA      2311      997407   1057     96264  24625
## 5          NA      2311      997407   1057     96264  24625
## 6          NA      2311      997407   1057     96264  24625
##   grassland_shrub   water wetlands agriculture_scaled high_kg_buff_scaled
## 1               0   13.44        0            -0.3525                  NA
## 2               0   13.44        0            -0.3525                  NA
## 3               0   13.44        0            -0.3525                  NA
## 4            8417 8353.56     1241             3.5449                  NA
## 5            8417 8353.56     1241             3.5449                  NA
## 6            8417 8353.56     1241             3.5449                  NA
names(BBS_ln) <- toupper(names(BBS_ln))
sort(unique(BBS_ln$DIET))
##  [1]               birds         carrion       fish          fruit        
##  [6] insects       mammals       nectar        omnivore      plants       
## [11] seeds         small_animals
## 12 Levels:  birds carrion fish fruit insects mammals nectar ... small_animals
BBS_ln$DIET <- relevel(BBS_ln$DIET, ref = "insects")
BBS_ln$YEAR_SCALED <- scale(BBS_ln$YEAR)
BBS_ln$FOREST_SCALED <- scale(BBS_ln$FOREST)
BBS_ln$AGRICULTURE_SCALED <- scale(BBS_ln$AGRICULTURE)
BBS_ln$FOREST_PROP <- BBS_ln$FOREST/BBS_ln$AGRICULTURE
BBS_ln$HIGH_KG_BUFF <- scale(BBS_ln$HIGH_KG_BUFF)
BBS_ln$DIET_VULN <- ifelse(BBS_ln$DIET %in% c("insects", "seeds", 'nectar', "birds"), "diet_vuln", "diet_invuln")
BBS_ln$DIET_VULN <- factor(BBS_ln$DIET_VULN)

BBS_ln$DIET_INSECTS <- ifelse(BBS_ln$DIET == "insects", "diet_insects", "diet_not_insects")
BBS_ln$DIET_INSECTS <- factor(BBS_ln$DIET_INSECTS)

with(BBS_ln, cor(FOREST, AGRICULTURE))
## [1] 0.262
# test_mod <- glm(SUM_ROUTE_ABUNDANCE ~ YEAR_SCALED + 
#                   AGRICULTURE_SCALED:HIGH_KG_BUFF_SCALED + 
#                   AGRICULTURE_SCALED +
#                   DIET_VULN:YEAR_SCALED + 
#                   DIET_VULN:HIGH_KG_BUFF_SCALED:YEAR_SCALED:AGRICULTURE_SCALED + 
#                   HIGH_KG_BUFF_SCALED:COMPOUND + 
#                   HIGH_KG_BUFF_SCALED, 
#                 data = BBS_ln[BBS_ln$COMPOUND %in% c("ATRAZINE", "GLYPHOSATE"),], 
#                 family = "poisson")
# summary(test_mod)
# #coefplot2(test_mod)
# #require(effects)
# 
# 
# test_mod2 <- glm(SUM_ROUTE_ABUNDANCE ~ YEAR_SCALED + 
#                   AGRICULTURE_SCALED:HIGH_KG_BUFF_SCALED + 
#                   AGRICULTURE_SCALED +
#                   DIET_INSECTS:YEAR_SCALED + 
#                   DIET_INSECTS:HIGH_KG_BUFF_SCALED:YEAR_SCALED:AGRICULTURE_SCALED + 
#                   HIGH_KG_BUFF_SCALED:COMPOUND + 
#                   HIGH_KG_BUFF_SCALED, 
#                 data = BBS_ln[BBS_ln$COMPOUND %in% c("ATRAZINE", "GLYPHOSATE"),], 
#                 family = "poisson")
# summary(test_mod2)
# #coefplot2(test_mod2)
# 
# 
# glmer_mod_3 <- glmer(SUM_ROUTE_ABUNDANCE ~ YEAR_SCALED + 
#                   AGRICULTURE_SCALED:HIGH_KG_BUFF_SCALED + 
#                   AGRICULTURE_SCALED +
#                   DIET_INSECTS:YEAR_SCALED + 
#                   DIET_INSECTS:HIGH_KG_BUFF_SCALED:YEAR_SCALED:AGRICULTURE_SCALED + 
#                   HIGH_KG_BUFF_SCALED:COMPOUND + 
#                   HIGH_KG_BUFF_SCALED +
#                     (1|ROUTE), 
#                 data = BBS_ln[BBS_ln$COMPOUND %in% c("ATRAZINE", "GLYPHOSATE"),], 
#                 family = "poisson")
# #coefplot2(glmer_mod_3)
# #coefplot2(list(glmer_mod_3, test_mod2))
# 
# glmer_mod_4 <- glmer(SUM_ROUTE_ABUNDANCE ~ YEAR_SCALED + 
#                   AGRICULTURE_SCALED:HIGH_KG_BUFF_SCALED + 
#                   AGRICULTURE_SCALED +
#                   DIET_INSECTS:YEAR_SCALED + 
#                     DIET_INSECTS +
#                   DIET_INSECTS:HIGH_KG_BUFF_SCALED:YEAR_SCALED:AGRICULTURE_SCALED + 
#                   HIGH_KG_BUFF_SCALED:COMPOUND + 
#                   HIGH_KG_BUFF_SCALED +
#                     (1|ROUTE), 
#                 data = BBS_ln[BBS_ln$COMPOUND %in% c("ATRAZINE", "GLYPHOSATE"),], 
#                 family = "poisson")
# summary(glmer_mod_4)
# #coefplot2(glmer_mod_4)
# #coefplot2(list(glmer_mod_3, test_mod2, glmer_mod_4))
# 
# glmer_mod_5 <- glmer(SUM_ROUTE_ABUNDANCE ~ 
#                        AGRICULTURE_SCALED:DIET_INSECTS + 
#                        AGRICULTURE_SCALED:HIGH_KG_BUFF_SCALED + 
#                        FOREST_SCALED:DIET_INSECTS +
#                        DIET_INSECTS:HIGH_KG_BUFF_SCALED:AGRICULTURE_SCALED + 
#                        HIGH_KG_BUFF_SCALED:COMPOUND:DIET_INSECTS + 
#                        HIGH_KG_BUFF_SCALED:DIET_INSECTS +
#                        (ROUTE|YEAR), 
#                 data = BBS_ln[BBS_ln$COMPOUND %in% c("ATRAZINE", "GLYPHOSATE"),], 
#                 family = "poisson")
# summary(glmer_mod_5)
#coefplot2(glmer_mod_5)
#coefplot2(list(glmer_mod_3, test_mod2, glmer_mod_4, glmer_mod_5))

Aggregating based on diet and then modelling

head(BBS_ln)
##   YEAR RTENO BUFFER  AOU ROUTEDATAID COUNTRYNUM STATENUM ROUTE RPID
## 1 1997 66001   1000 1320     6228264        840       66     1  101
## 2 1997 66001   1000 1320     6228264        840       66     1  101
## 3 1997 66001   1000 1320     6228264        840       66     1  101
## 4 1997 66001  10000 1320     6228264        840       66     1  101
## 5 1997 66001  10000 1320     6228264        840       66     1  101
## 6 1997 66001  10000 1320     6228264        840       66     1  101
##   SUM_ROUTE_ABUNDANCE    SCIENTIFIC_NAME X4_LETTER_CODE X6_LETTER_CODE
## 1                   2 ANAS PLATYRHYNCHOS           MALL         ANAPLA
## 2                   2 ANAS PLATYRHYNCHOS           MALL         ANAPLA
## 3                   2 ANAS PLATYRHYNCHOS           MALL         ANAPLA
## 4                   2 ANAS PLATYRHYNCHOS           MALL         ANAPLA
## 5                   2 ANAS PLATYRHYNCHOS           MALL         ANAPLA
## 6                   2 ANAS PLATYRHYNCHOS           MALL         ANAPLA
##     HABITAT  DIET NESTING BEHAVIOR CONSERVATION COMMON_NAME        ORDER
## 1 lake_pond seeds  ground  dabbler           LC     Mallard Anseriformes
## 2 lake_pond seeds  ground  dabbler           LC     Mallard Anseriformes
## 3 lake_pond seeds  ground  dabbler           LC     Mallard Anseriformes
## 4 lake_pond seeds  ground  dabbler           LC     Mallard Anseriformes
## 5 lake_pond seeds  ground  dabbler           LC     Mallard Anseriformes
## 6 lake_pond seeds  ground  dabbler           LC     Mallard Anseriformes
##     FAMILY GENUS       SPECIES     COMPOUND HIGH_KG_BUFF LOW_KG_BUFF
## 1 Anatidae  Anas platyrhynchos     ATRAZINE           NA          NA
## 2 Anatidae  Anas platyrhynchos   GLYPHOSATE           NA          NA
## 3 Anatidae  Anas platyrhynchos IMIDACLOPRID           NA          NA
## 4 Anatidae  Anas platyrhynchos     ATRAZINE           NA          NA
## 5 Anatidae  Anas platyrhynchos   GLYPHOSATE           NA          NA
## 6 Anatidae  Anas platyrhynchos IMIDACLOPRID           NA          NA
##   AG_PIX_BUFF BBS_ROUTE AGRICULTURE BARREN DEVELOPED FOREST
## 1          NA      2311       75241      0      7387      0
## 2          NA      2311       75241      0      7387      0
## 3          NA      2311       75241      0      7387      0
## 4          NA      2311      997407   1057     96264  24625
## 5          NA      2311      997407   1057     96264  24625
## 6          NA      2311      997407   1057     96264  24625
##   GRASSLAND_SHRUB   WATER WETLANDS AGRICULTURE_SCALED HIGH_KG_BUFF_SCALED
## 1               0   13.44        0            -0.3525                  NA
## 2               0   13.44        0            -0.3525                  NA
## 3               0   13.44        0            -0.3525                  NA
## 4            8417 8353.56     1241             3.5449                  NA
## 5            8417 8353.56     1241             3.5449                  NA
## 6            8417 8353.56     1241             3.5449                  NA
##   YEAR_SCALED FOREST_SCALED FOREST_PROP DIET_VULN     DIET_INSECTS
## 1      -1.882       -0.6396     0.00000 diet_vuln diet_not_insects
## 2      -1.882       -0.6396     0.00000 diet_vuln diet_not_insects
## 3      -1.882       -0.6396     0.00000 diet_vuln diet_not_insects
## 4      -1.882       -0.5070     0.02469 diet_vuln diet_not_insects
## 5      -1.882       -0.5070     0.02469 diet_vuln diet_not_insects
## 6      -1.882       -0.5070     0.02469 diet_vuln diet_not_insects
names(BBS_ln)
##  [1] "YEAR"                "RTENO"               "BUFFER"             
##  [4] "AOU"                 "ROUTEDATAID"         "COUNTRYNUM"         
##  [7] "STATENUM"            "ROUTE"               "RPID"               
## [10] "SUM_ROUTE_ABUNDANCE" "SCIENTIFIC_NAME"     "X4_LETTER_CODE"     
## [13] "X6_LETTER_CODE"      "HABITAT"             "DIET"               
## [16] "NESTING"             "BEHAVIOR"            "CONSERVATION"       
## [19] "COMMON_NAME"         "ORDER"               "FAMILY"             
## [22] "GENUS"               "SPECIES"             "COMPOUND"           
## [25] "HIGH_KG_BUFF"        "LOW_KG_BUFF"         "AG_PIX_BUFF"        
## [28] "BBS_ROUTE"           "AGRICULTURE"         "BARREN"             
## [31] "DEVELOPED"           "FOREST"              "GRASSLAND_SHRUB"    
## [34] "WATER"               "WETLANDS"            "AGRICULTURE_SCALED" 
## [37] "HIGH_KG_BUFF_SCALED" "YEAR_SCALED"         "FOREST_SCALED"      
## [40] "FOREST_PROP"         "DIET_VULN"           "DIET_INSECTS"
BBS_ln$HIGH_KG_BUFF <- as.numeric(BBS_ln$HIGH_KG_BUFF)
BBS_ln$AGRICULTURE_SCALED <- as.numeric(BBS_ln$AGRICULTURE_SCALED)
BBS_ln$HIGH_KG_BUFF_SCALED <- as.numeric(BBS_ln$HIGH_KG_BUFF_SCALED)
BBS_ln$YEAR_SCALED <- as.numeric(BBS_ln$YEAR_SCALED)
BBS_ln$FOREST_SCALED <- as.numeric(BBS_ln$FOREST_SCALED)
BBS_ln_diet <- data.frame(BBS_ln %>% group_by(YEAR, RTENO, BBS_ROUTE, BUFFER, ROUTE, COMPOUND, HIGH_KG_BUFF, AGRICULTURE, FOREST, DEVELOPED, AGRICULTURE_SCALED, YEAR_SCALED, FOREST_SCALED, HIGH_KG_BUFF_SCALED, FOREST_PROP, DIET) %>%
  summarise(SUM_ROUTE_ABUNDANCE_DIET = sum(SUM_ROUTE_ABUNDANCE)))
head(BBS_ln_diet)
##   YEAR RTENO BBS_ROUTE BUFFER ROUTE COMPOUND HIGH_KG_BUFF AGRICULTURE
## 1 1995 66061      2342    200    61 ATRAZINE      -0.2387       12628
## 2 1995 66061      2342    200    61 ATRAZINE      -0.2387       12628
## 3 1995 66061      2342    200    61 ATRAZINE      -0.2387       12628
## 4 1995 66061      2342    200    61 ATRAZINE      -0.2387       12628
## 5 1995 66061      2342    200    61 ATRAZINE      -0.2387       12628
## 6 1995 66061      2342    200    61 ATRAZINE      -0.2387       12628
##   FOREST DEVELOPED AGRICULTURE_SCALED YEAR_SCALED FOREST_SCALED
## 1   2330      1241            -0.6171       -2.37       -0.6271
## 2   2330      1241            -0.6171       -2.37       -0.6271
## 3   2330      1241            -0.6171       -2.37       -0.6271
## 4   2330      1241            -0.6171       -2.37       -0.6271
## 5   2330      1241            -0.6171       -2.37       -0.6271
## 6   2330      1241            -0.6171       -2.37       -0.6271
##   HIGH_KG_BUFF_SCALED FOREST_PROP     DIET SUM_ROUTE_ABUNDANCE_DIET
## 1             -0.2387      0.1845  insects                     1460
## 2             -0.2387      0.1845  carrion                        3
## 3             -0.2387      0.1845     fish                        3
## 4             -0.2387      0.1845    fruit                       40
## 5             -0.2387      0.1845   nectar                        1
## 6             -0.2387      0.1845 omnivore                      312
class(BBS_ln_diet$ROUTE)
## [1] "integer"
BBS_ln_diet$YEAR_fac <- as.factor(BBS_ln_diet$YEAR)
BBS_ln_diet$ROUTE <- as.factor(BBS_ln_diet$ROUTE)


range(na.omit(BBS_ln_diet$HIGH_KG_BUFF_SCALED))
## [1] -0.2858 14.3051
range(na.omit(BBS_ln_diet$AGRICULTURE_SCALED))
## [1] -0.6705  4.4323
range(na.omit(BBS_ln_diet$FOREST_SCALED))
## [1] -0.6396  4.7383
range(na.omit(BBS_ln_diet$FOREST_SCALED))
## [1] -0.6396  4.7383
class(BBS_ln_diet$DIET)
## [1] "factor"
class(BBS_ln_diet$ROUTE)
## [1] "factor"
range(BBS_ln_diet$SUM_ROUTE_ABUNDANCE_DIET)
## [1]    1 1500
ggplot(BBS_ln_diet, aes(x = ROUTE, y = SUM_ROUTE_ABUNDANCE_DIET)) +
  geom_point() +
  geom_boxplot() + theme_bw() +
  facet_wrap(~ YEAR)

plot of chunk unnamed-chunk-17

BBS_ln_diet$sqrt_response <- round(sqrt(BBS_ln_diet$SUM_ROUTE_ABUNDANCE_DIET), 0)
range(BBS_ln_diet$sqrt_response)
## [1]  1 39
BBS_ln_dietcc <- BBS_ln_diet[complete.cases(BBS_ln_diet),]
str(BBS_ln_dietcc)
## 'data.frame':    178487 obs. of  19 variables:
##  $ YEAR                    : int  1995 1995 1995 1995 1995 1995 1995 1995 1995 1995 ...
##  $ RTENO                   : int  66061 66061 66061 66061 66061 66061 66061 66061 66061 66061 ...
##  $ BBS_ROUTE               : int  2342 2342 2342 2342 2342 2342 2342 2342 2342 2342 ...
##  $ BUFFER                  : int  200 200 200 200 200 200 200 200 200 200 ...
##  $ ROUTE                   : Factor w/ 66 levels "1","2","3","5",..: 30 30 30 30 30 30 30 30 30 30 ...
##  $ COMPOUND                : Factor w/ 8 levels "ACETAMIPRID",..: 2 2 2 2 2 2 2 2 5 5 ...
##  $ HIGH_KG_BUFF            : num  -0.239 -0.239 -0.239 -0.239 -0.239 ...
##  $ AGRICULTURE             : num  12628 12628 12628 12628 12628 ...
##  $ FOREST                  : num  2330 2330 2330 2330 2330 ...
##  $ DEVELOPED               : num  1241 1241 1241 1241 1241 ...
##  $ AGRICULTURE_SCALED      : num  -0.617 -0.617 -0.617 -0.617 -0.617 ...
##  $ YEAR_SCALED             : num  -2.37 -2.37 -2.37 -2.37 -2.37 ...
##  $ FOREST_SCALED           : num  -0.627 -0.627 -0.627 -0.627 -0.627 ...
##  $ HIGH_KG_BUFF_SCALED     : num  -0.239 -0.239 -0.239 -0.239 -0.239 ...
##  $ FOREST_PROP             : num  0.184 0.184 0.184 0.184 0.184 ...
##  $ DIET                    : Factor w/ 12 levels "insects","","birds",..: 1 4 5 6 8 9 10 11 1 4 ...
##  $ SUM_ROUTE_ABUNDANCE_DIET: num  1460 3 3 40 1 312 22 337 1460 3 ...
##  $ YEAR_fac                : Factor w/ 17 levels "1995","1996",..: 1 1 1 1 1 1 1 1 1 1 ...
##  $ sqrt_response           : num  38 2 2 6 1 18 5 18 38 2 ...
load("predicting_glmers.Rdata")

# glmer_mod_6 <- glmer(na.action = na.omit,
#                      sqrt_response ~ 
#                        AGRICULTURE_SCALED:DIET+ 
#                        AGRICULTURE_SCALED:HIGH_KG_BUFF_SCALED + 
#                        FOREST_SCALED:DIET +
#                        DIET:HIGH_KG_BUFF_SCALED:AGRICULTURE_SCALED + 
#                        HIGH_KG_BUFF_SCALED:COMPOUND:DIET + 
#                        HIGH_KG_BUFF_SCALED:DIET +
#                        (1|YEAR_fac) + (1|ROUTE), 
#                      data = BBS_ln_dietcc[BBS_ln_dietcc$COMPOUND %in% c("ATRAZINE", 
#                                                                         "GLYPHOSATE")
#                                           & BBS_ln_dietcc$YEAR > 1996
#                                           & BBS_ln_dietcc$DIET != "",], 
#                      family = "poisson", 
#                      verbose = TRUE)
# sum_glmer_mod_6<-summary(glmer_mod_6)
# cf2_glmer_mod_6 <- coefplot2(glmer_mod_6)
# save(glmer_mod_6, sum_glmer_mod_6, cf2_glmer_mod_6, file = "glmer_mod_6.Rdata")
#load( "glmer_mod_6.Rdata")
#save.image("predicting_glmers.Rdata")
diet_levels <- levels(BBS_ln_dietcc$DIET)[-2]
ep <- expand.grid(AGRICULTURE_SCALED = seq(from = min(BBS_ln_dietcc$AGRICULTURE_SCALED), 
                                    to = max(BBS_ln_dietcc$AGRICULTURE_SCALED),
                                    by = 0.4),
                  FOREST_SCALED = mean(BBS_ln_dietcc$FOREST_SCALED),
                  DIET = diet_levels, 
                  HIGH_KG_BUFF_SCALED =seq(from = min(BBS_ln_dietcc$HIGH_KG_BUFF_SCALED), 
                                    to = max(BBS_ln_dietcc$HIGH_KG_BUFF_SCALED),
                                    by = 2),
                  COMPOUND = c("GLYPHOSATE", "ATRAZINE")
                    )

mm <- model.matrix(delete.response(terms(glmer_mod_6)), ep)

ep <- data.frame(ep, response = mm %*% fixef(glmer_mod_6))
ep <- with(ep, data.frame(ep, transformed.response = exp(response)))

pvar1 <- diag(mm %*% tcrossprod(vcov(glmer_mod_6), mm))
tvar1 <- pvar1 + VarCorr(glmer_mod_6)$YEAR + VarCorr(glmer_mod_6)$ROUTE
ep <- data.frame(ep, pse=sqrt(pvar1), tse=sqrt(tvar1))
ep <- with(ep,
                data.frame(ep,
                           plo=exp(response-1.96*pse),
                           phi=exp(response+1.96*pse),
                           tlo=exp(response-1.96*tse),
                           thi=exp(response+1.96*tse)))
head(ep)
##   AGRICULTURE_SCALED FOREST_SCALED    DIET HIGH_KG_BUFF_SCALED   COMPOUND
## 1            -0.6721     -0.004414 insects             -0.2818 GLYPHOSATE
## 2            -0.2721     -0.004414 insects             -0.2818 GLYPHOSATE
## 3             0.1279     -0.004414 insects             -0.2818 GLYPHOSATE
## 4             0.5279     -0.004414 insects             -0.2818 GLYPHOSATE
## 5             0.9279     -0.004414 insects             -0.2818 GLYPHOSATE
## 6             1.3279     -0.004414 insects             -0.2818 GLYPHOSATE
##   response transformed.response     pse    tse    plo    phi    tlo    thi
## 1    2.679               14.577 0.02996 0.1735 13.746 15.458 10.374 20.483
## 2    2.478               11.923 0.02981 0.1735 11.247 12.640  8.486 16.753
## 3    2.278                9.753 0.02981 0.1735  9.199 10.339  6.941 13.703
## 4    2.077                7.977 0.02998 0.1735  7.522  8.460  5.677 11.209
## 5    1.876                6.525 0.03030 0.1736  6.149  6.924  4.643  9.170
## 6    1.675                5.337 0.03077 0.1737  5.025  5.669  3.797  7.502
names(ep) <- tolower(names(ep))
names(ep) 
##  [1] "agriculture_scaled"   "forest_scaled"        "diet"                
##  [4] "high_kg_buff_scaled"  "compound"             "response"            
##  [7] "transformed.response" "pse"                  "tse"                 
## [10] "plo"                  "phi"                  "tlo"                 
## [13] "thi"
ep1 <- ep[ep$diet %in% c("insects", "omnivore"), ]
require(akima)
## Loading required package: akima
ep2 <- with(ep1[ep1$compound == "GLYPHOSATE" & ep1$diet == "insects",], interp(agriculture_scaled, high_kg_buff_scaled, transformed.response))
ep3 <- with(ep1[ep1$compound == "GLYPHOSATE" & ep1$diet == "omnivore",], interp(agriculture_scaled, high_kg_buff_scaled, transformed.response))

ep4 <- with(ep1[ep1$compound == "ATRAZINE" & ep1$diet == "insects",], interp(agriculture_scaled, high_kg_buff_scaled, transformed.response))
ep5 <- with(ep1[ep1$compound == "ATRAZINE" & ep1$diet == "omnivore",], interp(agriculture_scaled, high_kg_buff_scaled, transformed.response))


with(ep2, persp(x, y, z, theta = 330, phi = 20, 
                xlab = "Agriculture",
                ylab = "Pesticide application rate", 
                zlab = "Bird response",
                col = "royalblue",
                #nticks = 5, 
                #ticktype = "detailed",
                shade = 0.75))

plot of chunk unnamed-chunk-17

with(ep3, persp(x, y, z, theta = 55, phi = 30, 
                xlab = "Agriculture",
                ylab = "Pesticide application rate", 
                zlab = "Bird response",
                col = "cornflowerblue",
                #nticks = 5, 
                #ticktype = "detailed",
                shade = 0.75))

plot of chunk unnamed-chunk-17

with(ep4, persp(x, y, z, theta = 280, phi = 30, 
                xlab = "Agriculture",
                ylab = "Pesticide application rate", 
                zlab = "Bird response",
                col = "orange",
                #nticks = 5, 
                #ticktype = "detailed",
                shade = 0.75))

plot of chunk unnamed-chunk-17

with(ep5, persp(x, y, z, theta = 55, phi = 30, 
                xlab = "Agriculture",
                ylab = "Pesticide application rate", 
                zlab = "Bird response",
                col = "violetred",
                #nticks = 5, 
                #ticktype = "detailed",
                shade = 0.2))

plot of chunk unnamed-chunk-17

Earlier glm model which you can run if interested, heavily pseudo-replicated however.

glm_mod_6 <- glm(na.action = na.omit,
                 sqrt_response ~ 
                   DIET+
                   HIGH_KG_BUFF_SCALED +
                   AGRICULTURE_SCALED:DIET+ 
                   AGRICULTURE_SCALED:HIGH_KG_BUFF_SCALED + 
                   DIET:HIGH_KG_BUFF_SCALED:AGRICULTURE_SCALED + 
                   HIGH_KG_BUFF_SCALED:COMPOUND:DIET + 
                   HIGH_KG_BUFF_SCALED:DIET, 
                 data = BBS_ln_dietcc[BBS_ln_dietcc$COMPOUND %in% c("ATRAZINE", 
                                                                    "GLYPHOSATE")
                                      & BBS_ln_dietcc$YEAR > 1996
                                      & BBS_ln_dietcc$DIET != "",], 
                 family = "poisson")
summary(glm_mod_6)
## 
## Call:
## glm(formula = sqrt_response ~ DIET + HIGH_KG_BUFF_SCALED + AGRICULTURE_SCALED:DIET + 
##     AGRICULTURE_SCALED:HIGH_KG_BUFF_SCALED + DIET:HIGH_KG_BUFF_SCALED:AGRICULTURE_SCALED + 
##     HIGH_KG_BUFF_SCALED:COMPOUND:DIET + HIGH_KG_BUFF_SCALED:DIET, 
##     family = "poisson", data = BBS_ln_dietcc[BBS_ln_dietcc$COMPOUND %in% 
##         c("ATRAZINE", "GLYPHOSATE") & BBS_ln_dietcc$YEAR > 1996 & 
##         BBS_ln_dietcc$DIET != "", ], na.action = na.omit)
## 
## Deviance Residuals: 
##    Min      1Q  Median      3Q     Max  
## -3.124  -0.518  -0.143   0.452   4.097  
## 
## Coefficients:
##                                                           Estimate
## (Intercept)                                               3.175148
## DIETbirds                                                -3.134051
## DIETcarrion                                              -2.230065
## DIETfish                                                 -2.686719
## DIETfruit                                                -1.767611
## DIETmammals                                              -2.946772
## DIETnectar                                               -2.966412
## DIETomnivore                                             -0.731661
## DIETplants                                               -2.436557
## DIETseeds                                                -0.547610
## DIETsmall_animals                                        -2.846654
## HIGH_KG_BUFF_SCALED                                      -0.027365
## DIETinsects:AGRICULTURE_SCALED                            0.029820
## DIETbirds:AGRICULTURE_SCALED                              0.005821
## DIETcarrion:AGRICULTURE_SCALED                            0.084879
## DIETfish:AGRICULTURE_SCALED                               0.020911
## DIETfruit:AGRICULTURE_SCALED                              0.029860
## DIETmammals:AGRICULTURE_SCALED                            0.039498
## DIETnectar:AGRICULTURE_SCALED                            -0.017187
## DIETomnivore:AGRICULTURE_SCALED                           0.018079
## DIETplants:AGRICULTURE_SCALED                            -0.088285
## DIETseeds:AGRICULTURE_SCALED                              0.026454
## DIETsmall_animals:AGRICULTURE_SCALED                      0.016368
## HIGH_KG_BUFF_SCALED:AGRICULTURE_SCALED                    0.001428
## DIETbirds:HIGH_KG_BUFF_SCALED                             0.016325
## DIETcarrion:HIGH_KG_BUFF_SCALED                          -0.107665
## DIETfish:HIGH_KG_BUFF_SCALED                              0.030891
## DIETfruit:HIGH_KG_BUFF_SCALED                            -0.046663
## DIETmammals:HIGH_KG_BUFF_SCALED                          -0.030550
## DIETnectar:HIGH_KG_BUFF_SCALED                           -0.013445
## DIETomnivore:HIGH_KG_BUFF_SCALED                          0.006816
## DIETplants:HIGH_KG_BUFF_SCALED                            0.127757
## DIETseeds:HIGH_KG_BUFF_SCALED                             0.074079
## DIETsmall_animals:HIGH_KG_BUFF_SCALED                     0.004861
## DIETbirds:HIGH_KG_BUFF_SCALED:AGRICULTURE_SCALED         -0.000578
## DIETcarrion:HIGH_KG_BUFF_SCALED:AGRICULTURE_SCALED        0.010275
## DIETfish:HIGH_KG_BUFF_SCALED:AGRICULTURE_SCALED          -0.005944
## DIETfruit:HIGH_KG_BUFF_SCALED:AGRICULTURE_SCALED          0.002099
## DIETmammals:HIGH_KG_BUFF_SCALED:AGRICULTURE_SCALED        0.005676
## DIETnectar:HIGH_KG_BUFF_SCALED:AGRICULTURE_SCALED         0.006526
## DIETomnivore:HIGH_KG_BUFF_SCALED:AGRICULTURE_SCALED      -0.000291
## DIETplants:HIGH_KG_BUFF_SCALED:AGRICULTURE_SCALED         0.001589
## DIETseeds:HIGH_KG_BUFF_SCALED:AGRICULTURE_SCALED         -0.013408
## DIETsmall_animals:HIGH_KG_BUFF_SCALED:AGRICULTURE_SCALED -0.002292
## DIETinsects:HIGH_KG_BUFF_SCALED:COMPOUNDGLYPHOSATE        0.009158
## DIETbirds:HIGH_KG_BUFF_SCALED:COMPOUNDGLYPHOSATE          0.003118
## DIETcarrion:HIGH_KG_BUFF_SCALED:COMPOUNDGLYPHOSATE        0.034232
## DIETfish:HIGH_KG_BUFF_SCALED:COMPOUNDGLYPHOSATE          -0.000772
## DIETfruit:HIGH_KG_BUFF_SCALED:COMPOUNDGLYPHOSATE          0.021467
## DIETmammals:HIGH_KG_BUFF_SCALED:COMPOUNDGLYPHOSATE        0.007322
## DIETnectar:HIGH_KG_BUFF_SCALED:COMPOUNDGLYPHOSATE         0.002277
## DIETomnivore:HIGH_KG_BUFF_SCALED:COMPOUNDGLYPHOSATE       0.005580
## DIETplants:HIGH_KG_BUFF_SCALED:COMPOUNDGLYPHOSATE        -0.035643
## DIETseeds:HIGH_KG_BUFF_SCALED:COMPOUNDGLYPHOSATE         -0.005885
## DIETsmall_animals:HIGH_KG_BUFF_SCALED:COMPOUNDGLYPHOSATE  0.010567
##                                                          Std. Error
## (Intercept)                                                0.002811
## DIETbirds                                                  0.029237
## DIETcarrion                                                0.010397
## DIETfish                                                   0.012308
## DIETfruit                                                  0.008262
## DIETmammals                                                0.025355
## DIETnectar                                                 0.019745
## DIETomnivore                                               0.004942
## DIETplants                                                 0.017813
## DIETseeds                                                  0.004607
## DIETsmall_animals                                          0.014476
## HIGH_KG_BUFF_SCALED                                        0.005701
## DIETinsects:AGRICULTURE_SCALED                             0.005173
## DIETbirds:AGRICULTURE_SCALED                               0.053022
## DIETcarrion:AGRICULTURE_SCALED                             0.018927
## DIETfish:AGRICULTURE_SCALED                                0.021730
## DIETfruit:AGRICULTURE_SCALED                               0.014734
## DIETmammals:AGRICULTURE_SCALED                             0.049625
## DIETnectar:AGRICULTURE_SCALED                              0.039123
## DIETomnivore:AGRICULTURE_SCALED                            0.007512
## DIETplants:AGRICULTURE_SCALED                              0.030605
## DIETseeds:AGRICULTURE_SCALED                               0.006612
## DIETsmall_animals:AGRICULTURE_SCALED                       0.026318
## HIGH_KG_BUFF_SCALED:AGRICULTURE_SCALED                     0.001109
## DIETbirds:HIGH_KG_BUFF_SCALED                              0.059939
## DIETcarrion:HIGH_KG_BUFF_SCALED                            0.022122
## DIETfish:HIGH_KG_BUFF_SCALED                               0.024384
## DIETfruit:HIGH_KG_BUFF_SCALED                              0.017615
## DIETmammals:HIGH_KG_BUFF_SCALED                            0.057678
## DIETnectar:HIGH_KG_BUFF_SCALED                             0.045433
## DIETomnivore:HIGH_KG_BUFF_SCALED                           0.010036
## DIETplants:HIGH_KG_BUFF_SCALED                             0.031491
## DIETseeds:HIGH_KG_BUFF_SCALED                              0.009155
## DIETsmall_animals:HIGH_KG_BUFF_SCALED                      0.029825
## DIETbirds:HIGH_KG_BUFF_SCALED:AGRICULTURE_SCALED           0.011957
## DIETcarrion:HIGH_KG_BUFF_SCALED:AGRICULTURE_SCALED         0.004264
## DIETfish:HIGH_KG_BUFF_SCALED:AGRICULTURE_SCALED            0.004806
## DIETfruit:HIGH_KG_BUFF_SCALED:AGRICULTURE_SCALED           0.003469
## DIETmammals:HIGH_KG_BUFF_SCALED:AGRICULTURE_SCALED         0.011524
## DIETnectar:HIGH_KG_BUFF_SCALED:AGRICULTURE_SCALED          0.009166
## DIETomnivore:HIGH_KG_BUFF_SCALED:AGRICULTURE_SCALED        0.001951
## DIETplants:HIGH_KG_BUFF_SCALED:AGRICULTURE_SCALED          0.005777
## DIETseeds:HIGH_KG_BUFF_SCALED:AGRICULTURE_SCALED           0.001786
## DIETsmall_animals:HIGH_KG_BUFF_SCALED:AGRICULTURE_SCALED   0.005722
## DIETinsects:HIGH_KG_BUFF_SCALED:COMPOUNDGLYPHOSATE         0.003159
## DIETbirds:HIGH_KG_BUFF_SCALED:COMPOUNDGLYPHOSATE           0.032824
## DIETcarrion:HIGH_KG_BUFF_SCALED:COMPOUNDGLYPHOSATE         0.012106
## DIETfish:HIGH_KG_BUFF_SCALED:COMPOUNDGLYPHOSATE            0.013002
## DIETfruit:HIGH_KG_BUFF_SCALED:COMPOUNDGLYPHOSATE           0.009379
## DIETmammals:HIGH_KG_BUFF_SCALED:COMPOUNDGLYPHOSATE         0.033707
## DIETnectar:HIGH_KG_BUFF_SCALED:COMPOUNDGLYPHOSATE          0.025555
## DIETomnivore:HIGH_KG_BUFF_SCALED:COMPOUNDGLYPHOSATE        0.004578
## DIETplants:HIGH_KG_BUFF_SCALED:COMPOUNDGLYPHOSATE          0.014605
## DIETseeds:HIGH_KG_BUFF_SCALED:COMPOUNDGLYPHOSATE           0.003939
## DIETsmall_animals:HIGH_KG_BUFF_SCALED:COMPOUNDGLYPHOSATE   0.016377
##                                                          z value Pr(>|z|)
## (Intercept)                                              1129.49  < 2e-16
## DIETbirds                                                -107.20  < 2e-16
## DIETcarrion                                              -214.49  < 2e-16
## DIETfish                                                 -218.29  < 2e-16
## DIETfruit                                                -213.95  < 2e-16
## DIETmammals                                              -116.22  < 2e-16
## DIETnectar                                               -150.24  < 2e-16
## DIETomnivore                                             -148.06  < 2e-16
## DIETplants                                               -136.79  < 2e-16
## DIETseeds                                                -118.86  < 2e-16
## DIETsmall_animals                                        -196.65  < 2e-16
## HIGH_KG_BUFF_SCALED                                        -4.80  1.6e-06
## DIETinsects:AGRICULTURE_SCALED                              5.76  8.2e-09
## DIETbirds:AGRICULTURE_SCALED                                0.11   0.9126
## DIETcarrion:AGRICULTURE_SCALED                              4.48  7.3e-06
## DIETfish:AGRICULTURE_SCALED                                 0.96   0.3359
## DIETfruit:AGRICULTURE_SCALED                                2.03   0.0427
## DIETmammals:AGRICULTURE_SCALED                              0.80   0.4261
## DIETnectar:AGRICULTURE_SCALED                              -0.44   0.6604
## DIETomnivore:AGRICULTURE_SCALED                             2.41   0.0161
## DIETplants:AGRICULTURE_SCALED                              -2.88   0.0039
## DIETseeds:AGRICULTURE_SCALED                                4.00  6.3e-05
## DIETsmall_animals:AGRICULTURE_SCALED                        0.62   0.5340
## HIGH_KG_BUFF_SCALED:AGRICULTURE_SCALED                      1.29   0.1978
## DIETbirds:HIGH_KG_BUFF_SCALED                               0.27   0.7853
## DIETcarrion:HIGH_KG_BUFF_SCALED                            -4.87  1.1e-06
## DIETfish:HIGH_KG_BUFF_SCALED                                1.27   0.2052
## DIETfruit:HIGH_KG_BUFF_SCALED                              -2.65   0.0081
## DIETmammals:HIGH_KG_BUFF_SCALED                            -0.53   0.5963
## DIETnectar:HIGH_KG_BUFF_SCALED                             -0.30   0.7673
## DIETomnivore:HIGH_KG_BUFF_SCALED                            0.68   0.4970
## DIETplants:HIGH_KG_BUFF_SCALED                              4.06  5.0e-05
## DIETseeds:HIGH_KG_BUFF_SCALED                               8.09  5.9e-16
## DIETsmall_animals:HIGH_KG_BUFF_SCALED                       0.16   0.8705
## DIETbirds:HIGH_KG_BUFF_SCALED:AGRICULTURE_SCALED           -0.05   0.9614
## DIETcarrion:HIGH_KG_BUFF_SCALED:AGRICULTURE_SCALED          2.41   0.0159
## DIETfish:HIGH_KG_BUFF_SCALED:AGRICULTURE_SCALED            -1.24   0.2162
## DIETfruit:HIGH_KG_BUFF_SCALED:AGRICULTURE_SCALED            0.61   0.5451
## DIETmammals:HIGH_KG_BUFF_SCALED:AGRICULTURE_SCALED          0.49   0.6223
## DIETnectar:HIGH_KG_BUFF_SCALED:AGRICULTURE_SCALED           0.71   0.4764
## DIETomnivore:HIGH_KG_BUFF_SCALED:AGRICULTURE_SCALED        -0.15   0.8816
## DIETplants:HIGH_KG_BUFF_SCALED:AGRICULTURE_SCALED           0.28   0.7832
## DIETseeds:HIGH_KG_BUFF_SCALED:AGRICULTURE_SCALED           -7.51  6.1e-14
## DIETsmall_animals:HIGH_KG_BUFF_SCALED:AGRICULTURE_SCALED   -0.40   0.6888
## DIETinsects:HIGH_KG_BUFF_SCALED:COMPOUNDGLYPHOSATE          2.90   0.0037
## DIETbirds:HIGH_KG_BUFF_SCALED:COMPOUNDGLYPHOSATE            0.09   0.9243
## DIETcarrion:HIGH_KG_BUFF_SCALED:COMPOUNDGLYPHOSATE          2.83   0.0047
## DIETfish:HIGH_KG_BUFF_SCALED:COMPOUNDGLYPHOSATE            -0.06   0.9526
## DIETfruit:HIGH_KG_BUFF_SCALED:COMPOUNDGLYPHOSATE            2.29   0.0221
## DIETmammals:HIGH_KG_BUFF_SCALED:COMPOUNDGLYPHOSATE          0.22   0.8280
## DIETnectar:HIGH_KG_BUFF_SCALED:COMPOUNDGLYPHOSATE           0.09   0.9290
## DIETomnivore:HIGH_KG_BUFF_SCALED:COMPOUNDGLYPHOSATE         1.22   0.2229
## DIETplants:HIGH_KG_BUFF_SCALED:COMPOUNDGLYPHOSATE          -2.44   0.0147
## DIETseeds:HIGH_KG_BUFF_SCALED:COMPOUNDGLYPHOSATE           -1.49   0.1352
## DIETsmall_animals:HIGH_KG_BUFF_SCALED:COMPOUNDGLYPHOSATE    0.65   0.5188
##                                                             
## (Intercept)                                              ***
## DIETbirds                                                ***
## DIETcarrion                                              ***
## DIETfish                                                 ***
## DIETfruit                                                ***
## DIETmammals                                              ***
## DIETnectar                                               ***
## DIETomnivore                                             ***
## DIETplants                                               ***
## DIETseeds                                                ***
## DIETsmall_animals                                        ***
## HIGH_KG_BUFF_SCALED                                      ***
## DIETinsects:AGRICULTURE_SCALED                           ***
## DIETbirds:AGRICULTURE_SCALED                                
## DIETcarrion:AGRICULTURE_SCALED                           ***
## DIETfish:AGRICULTURE_SCALED                                 
## DIETfruit:AGRICULTURE_SCALED                             *  
## DIETmammals:AGRICULTURE_SCALED                              
## DIETnectar:AGRICULTURE_SCALED                               
## DIETomnivore:AGRICULTURE_SCALED                          *  
## DIETplants:AGRICULTURE_SCALED                            ** 
## DIETseeds:AGRICULTURE_SCALED                             ***
## DIETsmall_animals:AGRICULTURE_SCALED                        
## HIGH_KG_BUFF_SCALED:AGRICULTURE_SCALED                      
## DIETbirds:HIGH_KG_BUFF_SCALED                               
## DIETcarrion:HIGH_KG_BUFF_SCALED                          ***
## DIETfish:HIGH_KG_BUFF_SCALED                                
## DIETfruit:HIGH_KG_BUFF_SCALED                            ** 
## DIETmammals:HIGH_KG_BUFF_SCALED                             
## DIETnectar:HIGH_KG_BUFF_SCALED                              
## DIETomnivore:HIGH_KG_BUFF_SCALED                            
## DIETplants:HIGH_KG_BUFF_SCALED                           ***
## DIETseeds:HIGH_KG_BUFF_SCALED                            ***
## DIETsmall_animals:HIGH_KG_BUFF_SCALED                       
## DIETbirds:HIGH_KG_BUFF_SCALED:AGRICULTURE_SCALED            
## DIETcarrion:HIGH_KG_BUFF_SCALED:AGRICULTURE_SCALED       *  
## DIETfish:HIGH_KG_BUFF_SCALED:AGRICULTURE_SCALED             
## DIETfruit:HIGH_KG_BUFF_SCALED:AGRICULTURE_SCALED            
## DIETmammals:HIGH_KG_BUFF_SCALED:AGRICULTURE_SCALED          
## DIETnectar:HIGH_KG_BUFF_SCALED:AGRICULTURE_SCALED           
## DIETomnivore:HIGH_KG_BUFF_SCALED:AGRICULTURE_SCALED         
## DIETplants:HIGH_KG_BUFF_SCALED:AGRICULTURE_SCALED           
## DIETseeds:HIGH_KG_BUFF_SCALED:AGRICULTURE_SCALED         ***
## DIETsmall_animals:HIGH_KG_BUFF_SCALED:AGRICULTURE_SCALED    
## DIETinsects:HIGH_KG_BUFF_SCALED:COMPOUNDGLYPHOSATE       ** 
## DIETbirds:HIGH_KG_BUFF_SCALED:COMPOUNDGLYPHOSATE            
## DIETcarrion:HIGH_KG_BUFF_SCALED:COMPOUNDGLYPHOSATE       ** 
## DIETfish:HIGH_KG_BUFF_SCALED:COMPOUNDGLYPHOSATE             
## DIETfruit:HIGH_KG_BUFF_SCALED:COMPOUNDGLYPHOSATE         *  
## DIETmammals:HIGH_KG_BUFF_SCALED:COMPOUNDGLYPHOSATE          
## DIETnectar:HIGH_KG_BUFF_SCALED:COMPOUNDGLYPHOSATE           
## DIETomnivore:HIGH_KG_BUFF_SCALED:COMPOUNDGLYPHOSATE         
## DIETplants:HIGH_KG_BUFF_SCALED:COMPOUNDGLYPHOSATE        *  
## DIETseeds:HIGH_KG_BUFF_SCALED:COMPOUNDGLYPHOSATE            
## DIETsmall_animals:HIGH_KG_BUFF_SCALED:COMPOUNDGLYPHOSATE    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 501735  on 62855  degrees of freedom
## Residual deviance:  43072  on 62801  degrees of freedom
## AIC: 254471
## 
## Number of Fisher Scoring iterations: 4
coefplot2(glm_mod_6)

plot of chunk unnamed-chunk-18